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ABSTRACT 


The  ability  of  a  submarine  to  maintain  ordered  depth,  especially  during  periscope 
depth  operations  at  low  speeds,  is  vital  for  the  vessel  to  perform  its  mission  and  avoid 
detection  Modem  submarines  exhibit  an  inherent  phenomenon  that  produces  an 
undesirable  ship  response  at  low  speeds,  commonly  referred  to  as  dive  plane  reversal  The 
physical  parameters  that  govern  this  occurrence  are  related  in  this  thesis  to  the  problem  of 
multiple  steady  state  solutions  in  the  vertical  plane 

Generic  solution  branching,  in  the  form  of  pitchfork  bifurcations,  can  occur  when  the 
nominal  level  flight  path  loses  its  stability  A  systematic  study  reveals  the  existence  of  a 
critical  Froude  number,  based  on  the  vessel's  speed  and  metacentric  height,  where  this 
branching  occurs  Bifurcation  theory  techniques  and  numerical  computations  are  utilized 
to  classify  the  effect  that  geometric  parameters,  trim  and  ballast  conditions,  and 
hydrodynamic  properties  have  on  the  existence  of  these  multiple  solutions 
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I.  INTRODUCTION 


One  of  the  most  critical  functions  of  a  submersible  is  accurate  depth  keeping  at  the 
commanded  depth  Such  a  function  can  be  carried  out  either  manually  or  automatically, 
especially  in  cases  where  human  intervention  is  not  possible  or  desirable  Due  to  the 
technological  significance  and  numerous  scientific  applications  of  submersible  vehicle 
systems,  design  issues  of  appropriate  depth  keeping  control  laws  have  received  wide 
attention  in  the  past  Such  control  system  designs  include  linear  and  nonlinear  controllers 
[Refs  1&2],  model  based  compensators  [Ref  3],  adaptive  control  [Ref  4],  and  sliding 
mode  control  laws  [Refs  5&6] 

Response  accuracy  and  stability  are  primary  considerations  in  designing  depth 
keeping  control  law  Of  paramount  importance,  in  this  area,  are  the  robustness  properties 
of  the  particular  design,  i  e ,  its  ability  to  maintain  accuracy  and  stability  in  the  presence  of 
incomplete  sensor  and  environmental  information,  as  well  as  actual/mathematical  model 
mismatch  The  scope  of  the  work  in  this  thesis  is  to  demonstrate  a  potential  loss  of 
stability  that  may  occur  when  a  submarine  is  operating  at  low  speeds  This  loss  of  stability 
can  occur  regardless  of  the  particular  means  used  for  depth  control  The  study  is 
accomplished  through  the  use  of  an  eigenvalue  analysis,  a  steady  state  analysis  and  a 
controllability  analysis  [Ref  7],  It  is  shown  that  such  a  loss  of  stability  is  accompanied  by  a 
slow  divergence  of  trajectories  away  from  the  commanded  path  Solution  branching 
occurs  in  the  form  of  generic  pitchfork  bifurcations  [Refs  8,  9]  A  complete 
characterization  of  the  problem  is  given  utilizing  singularity  techniques,  which  have  been 
proven  to  be  very  useful  in  the  analysis  of  similar  problems  [Refs  10,  11,  12]  The  use  of 


bifurcation  theory  allows  the  crucial  vehicle  parameters  that  govern  the  problem  of 
solution  branching  to  be  determined,  and  the  develop  guidelines  to  prevent  its  occurrence 
Finally,  a  new  look  at  the  problem  of  dive  plane  reversal  [Ref  13],  based  on  solution 
branching  results  is  presented  The  term  dive  plane  reversal  refers  to  a  well-known 
phenomenon  in  submarine  operations  where,  .^ring  low  speed  depth  keeping,  there  is  a 
need  to  reverse  the  direction  of  dive  plane  defl*  ction  in  order  to  execute  a  given  change  in 
depth  Physically,  this  can  be  explained  by  considering  the  relative  magnitude  of  the 
hydrodynamic  forces  At  moderate  and  high  speeds,  the  normal  force  on  the  submarine's 
hull  due  to  the  angle  of  attack  exceeds  the  normal  I'ive  plane  force  and  the  boat  responds 
to  ordered  dive  plane  angles  as  expected.  The  phenomenon  of  dive  plane  reversal  occurs 
at  speeds  below  a  certain  critical  speed  in  which  the  normal  hull  force  is  less  than  the 
normal  dive  plane  force  and  the  response  of  the  boat  is  reversed.  [Ref  14]  Vehicle 
modeling  in  this  work  follows  standard  notation  [Ref  15]  and  numerical  results  are 
presented  for  the  DARPA  SUBOFF  model  [Ref  16]  for  which  a  set  of  hydrodynamic 
coefficients  and  geometric  properties  is  available  Special  emphasis  is  given  to  identifying 
the  proper  non-dimensional  parameters  in  the  problem,  so  that  extension  of  these  results 
to  full  scale  models  and  other  designs  is  possible  using  minimal  experimental  and/or 
analytical  results 


n.  VEHICLE  MODELING 


A.  EQUATIONS  OF  MOTION 

1.  Introduction 

For  the  purpose  of  this  problem,  and  the  subject  of  the  maneuvering  and  motion 
control  of  the  vehicle,  the  following  assumptions  are  made 

1 )  The  vehicle  behaves  as  a  rigid  body, 

2)  The  earth's  rotation  is  negligible  as  far  as  acceleration  components  of  the 
center  of  mass  are  concerned, 

3)  The  primary  forces  that  act  on  the  vehicle  have  inertial,  gravitational, 
hydrostatic  and  hydrodynamic  origins 

2.  Coordinate  Systems  And  Positional  Definitions 

A  global  navigation  frame,  OA77,  as  shown  in  Figure  2  1,  is  defined  with  origin, 
O,  and  a  set  of  axes  aligned  with  directions  North,  East  and  Down  This  produces  a  right- 
hand  reference  frame  with  unit  vectors  /.  J,  and  K.  Ignoring  the  earth's  rotation  rate  in 
comparison  to  the  angular  rates  produced  by  the  vessel's  motion,  it  can  be  said  that  the  7. 
J.  and  K  coordinate  frame  is  an  inertial  reference  frame  in  which  Newton's  Laws  of 
Motion  will  be  valid  As  seen  in  Figure  2  1,  a  vehicle's  position  ,  in  this  frame  will  have 
the  vector  components,  I  ^  standard  convention  that  is  used 

in  both  aircraft  and  marine  vehicle  dynamics  places  the  Y  axis  to  the  right  while  looking 
along  the  X  axis,  and  the  Z  axis  is  positive  downwards  Figure  2  1  also  shows  a  vehicle 
with  some  general  attitude  in  the  original  coordinate  system 


Now  define  a  body  fixed  frame  of  reference  Oryx,  with  the  origin  O ,  located  on 
the  vehicle  centerline,  moving  and  rotating  with  the  vehicle,  in  which  the  vehicle's  center 
of  mass,  G,  has  some  position  other  than  the  origin  of  the  vehicle  fixed  frame  (refer  to 
Figure  2-1)  This  origin  will  be  the  point  about  which  all  vehicle  body  force  will  be 
computed  in  later  sections  of  this  chapter  The  convention  used  for  the  vehicle  fixed 
frame,  with  unit  vectors  T ,  J ,  ^ ,  is  that  its  origin  lies  at  the  ship's  center  (origin  at  the  half 
length),  the  x  axis  at  main  deck  level,  parallel  to  the  longitudinal  centerline,  with  the  z  axis 
vertically  down  The  vessel's  center  of  gravity  and  center  of  buoyancy  do  not  generally  lie 
at  the  origin  of  the  body  fixed  frame,  nor  are  they  collocated  These  points  are  denoted  by 
G  and  B  respectively  The  vector  components  of  and  Pg  are  thus  [ xj  +  }'a]  /. 

and  f  xJ  +  }gj Zgic  J  The  location  of  the  center  of  mass  is  important  because 
Newton's  Laws  of  Motion  equate  forces  on  a  body  to  the  rate  of  change  of  linear 
momentum  of  the  center  of  mass  and  moments  about  the  body  center  of  mass  to  the  rate 
of  change  of  angular  momentum  This  is  a  particularly  important  point  to  bear  in  mind  for 
ship  and  submarine  motion  dynamics  as  the  center  of  buoyancy  is  determined  by  the  shape 
of  the  submerged  portion  of  the  ship  body  while  the  center  of  gravity  is  determined  by  the 
distribution  of  the  weight  over  the  entire  ship  body  Having  defined  a  coordinate  system 
that  will  be  used  to  describe  a  ship's  position,  there  is  a  need  to  define  angular  orientation 
of  the  ship's  body,  leading  to  the  definition  of  translational  and  rotational  velocities  and 
accelerations 

3.  Angular  Position  In  The  Global  Reference  Frame 

There  are  several  different  ways  that  the  attitude  of  a  vehicle  can  be  described  in 
reference  to  the  global  frame  The  most  usual  and  common  method  is  to  define  three 
angles,  called  Euler  angles  that  uniquely  define  the  angular  orientation  of  the  vehicle 
reference  frame,  relative  to  the  global  reference  frame  One  problem  with  the  Euler  angle 


•  • 


approach  is  that  a  singularity  exists  when  one  of  the  angles  reaches  90  degrees  (an  aircraft 
in  pure  vertical  flight  cannot  distinguish  its  azimuth  angle  from  its  roll  angle).  This 
limitation  which  can  sometimes,  although  rarely,  cause  trouble  in  flight  simulations  and 
control  computations  can  be  overcome  by  the  use  of  quaternions  which  introduce  four 
rather  than  three  variables  to  describe  angular  position.  In  this  presentation,  however,  we 
will  use  Euler  angles  as  it  is  the  most  widely  used  method,  although  the  use  of  quaternions 
has  found  favor  in  robotics  applications  and  computer  graphics. 

While  any  consistent  definition  of  three  base  angles  would  be  sufficient,  the  most 
convenient  formulation  of  these  angles  is  in  the  widely  used  military  terms  of  an  azimuth, 
elevation,  and  spin  notation  The  rates  of  change  of  these  angles  do  not  generally 
correspond  to  the  other  commonly  used  angular  rates  describing  angular  velocity  of  a 
vehicle,  that  is  yaw  rate,  pitch  rate  and  roll  rate  except,  as  will  be  seen,  where  motion  is 
limited  to  small  angle  rotations 

4.  Rotational  Transformations 

Vehicle  attitude  is  important  in  defining  position.  Under  dynamic  conditions,  the 
pitch  or  roll  can  cause  problems  for  manned  vessels  and  the  heading  is  critical  in 
navigation  For  the  purpose  of  considering  angular  rotations,  consider  a  forward 
transformation  from  a  coordinate  triad  aligned  with  the  global  reference  frame  and  perform 
a  sequence  of  three  rotations  to  finally  align  the  result  with  a  frame  that  is  assumed  to  be 
parallel  to  the  vehicle  body  coordinate  axes,  and  moving  with  the  vehicle  at  all  times 
Begin  by  defining  an  azimuth  rotation,  v ,  as  a  positive  rotation  about  the  global  Z  axis 
Next  define  a  subsequent  rotation  6 ,  (positive  up)  about  the  new  Y  axis,  followed  by  a 
positive  rotation  4i,  about  the  new  X  axis  The  triple  rotational  transformation  in  terms  of 
these  three  angles,  is  then  sufficient  to  describe  the  angular  orientation  of  the  vehicle  at  any 
time  It  follows  that  any  position  vector,  R^,  in  an  original  reference  frame  given  by 


K-  l^o-yo’^ol '  different  coordinates  in  a  rotated  frame  when  an  azimuth 

rotation  by  angle  V|t ,  is  made  about  the  global  Z  axis. 

If  the  new  position  is  defined  by,  ^  it  can  be  seen  that  there  is  a 

relation  between  the  vector's  coordinates  in  the  new  reference  frame  with  those  that  it  had 
in  the  old  reference  frame  Using  trigonometrical  relationships  and  Figure  2  2  as  a  guide,  it 
follows  that, 

X,  =  X^cos\/^Y^sm\if  (2  1) 

Y^=-X^  sin  y  +  };  cosy  (2  2) 

This  relation  can  be  expressed  in  matrix  form  by  the  rotation  matrix  operation, 

(2.3) 

where  the  rotation  matrix  [T^z]'  represents  an  orthogonal  transformation 
Premultiplication  of  this  rotation  matrix  with  any  vector,  will  produce  the  components 
of  the  same  vector  in  the  rotated  coordinate  frame.  Continuing  with  the  series  of  rotations 
results  in  a  combined  total  rotation  transformation, 

rr(t).0.y>i=TOTO7'ry;  (2  4) 

In  expanded  notation  equation  2.4  takes  the  form; 

cosy  cos 0  s/>iycos0  -sinQ 

cosy  s/«0s/n<|» -s/ny  cos <(>  s/ny  s/«0s/« (j)  + cosy  cos (])  cosQsin^  , 
cosy S//I0 cos ((n-s/nys/>i<(»  s/nys/w0cos(|>-cosys//?(l)  cos0cos<t> 

and  it  can  be  said  that  any  position  vector  in  an  original  reference  frame  may  be 

expressed  in  a  rotated  frame  with  coordinates  given  by  the  operation, 

^=[rr(t>.0.y>R«  (2  5) 


*j 


4r 
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Figure  2*2.  Azimuth  Rotation. 


5.  Kinematics 

Kinematics  deals  with  the  relationships  of  motion  quantities  regardless  of  the 
forces  induced  by  their  prescribed  motions  Descriptions  of  a  ship's  positiqn  both 
translational  and  rotational,  will  need  to  be  related  to  velocities,  both  translational  and 
rotational,  prior  to  e?aending  the  situation  to  accelerations  The  connection  between 
translational  velocity  and  the  rate  of  change  of  translational  position  is  straight  forward 
Define  the  global  velocity  as. 


h  = 


X 

Y 

Z 


(26) 


This  vector  will  have  components  that  are  different  when  seen  in  a  body  fixed  frame 
Define  the  body  fixed  components  of  the  global  velocity  vector  as  [u,  v,  w]'  These 
components,  in  terms  of  the  above  global  quantities  are  given  by  the  forward 
transformation  defined  earlier  to  be. 


u 

X 

V 

w 

L  J 

Y 

Z 

(2.7) 


It  is  a  simple  reverse  coordinate  transformation  from  body  fixed  to  global  coordinates  to 
see  that. 


X 

u 

Y 

V 

Z 

H 

This  shows  that  the  progression  of  a  vehicle  in  a  global  frame  clearly  depends  on  its  local 
velocity  components  and  its  attitude  Put  simply,  u  corresponding  to  a  vehicle's  forward 
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speed  (surge),  v  corresponding  to  a  sideslip  velocity  (sway),  and  w  corresponding  to  any 
velocity  component  in  the  local  Z  direction  (heave)  and  the  vehicle's  global  velocity 
components  depend  on  heading,  pitch,  and  roll  attitude 

The  connection  between  angular  attitude  and  angular  velocity  is  not  quite  so 
apparent  At  first  sight,  it  is  tempting  to  define  the  instantaneous  angular  velocity  of  the 
vehicle  simply  as  the  rate  of  change  of  its  angular  position  defined  by  the  Euler  angles 
This  is  erroneous  however,  because  the  rotation  8 ,  was  defined  as  a  rotation  about  the 
intermediate  frame  after  a  rotation  y  had  been  made  Vehicle  inertial  angular  rates  are 
defined  in  terms  of  components  that  have  angular  velocities  about  the  global  axes  It  is 
necessary  to  relate  both  Euler  angles  and  their  rates  of  change  to  angular  velocity 
components  about  the  global  axes  to  their  components  lying  along  the  body  fixed  axes  in 
any  attitude  The  prime  reason  for  this  is  that  it  is  difficult  to  construct  physical  sensors  to 
measure  rates  of  change  of  Euler  angles  However,  rate  gyros  in  common  use  today  are 
easily  constructed  to  measure  the  components  of  the  inertial  angular  velocity  of  a  vehicle 
that  lie  along  the  vehicle's  body  axes  It  follows  that  the  instantaneous  angular  velocity  of 
the  vehicle  can  be  related  to  the  instantaneous  rate  of  change  of  angular  orientation  only 
after  considerations  of  the  intermediate  transformations  used  In  other  words,  if  a  is 
defined  as  the  angular  attitude  vector,  a  =  [v,0,0],  and  the  inertial  components  of  the 
vehicle  angular  rate  lying  along  the  body  axes  as  to  -[p.q.r] ,  then,  d  =  The  details 

of  the  nonlinear  functional  relations  involved  are  provided  by  viewing  the  rate  of  change  of 
the  rotation  \|/  as  a  vector  quantity  lying  along  the  original  Z  axis  The  rate  of  change  of 
the  angle  8  is  viewed  as  a  vector  quantity  lying  along  the  Y  axis  of  the  first  intermediate 
frame,  and  the  rate  of  change  of  the  angle  d  is  viewed  as  a  vector  lying  along  the  X  axis  of 
the  finql  (body  fixed)  frame  Each  of  these  component  rates  of  change  of  angular  position 
has  component  parts  that  project  onto  the  final  frame  and  it  is  the  sum  total  of  all  the 
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components  that  give  the  total  angular  velocity  as  seen  in  the  final  frame  of  reference 
Using  the  required  transformations  for  the  rate  components  from  each  Euler  angle,  we  get. 


P 

<1 

r 


0 

+  e 
0 


+  0 
0 


with  the  result. 


(29) 


= 

vir5//j6+6coi((i 

r 

y  cosQcos  <j>  -  6  5//1 0 

Inverting  equation  2  10  yields  a  solution  for  the  rates  of  change  of  the  Euler  angles  in 
terms  of  the  body  fixed  components  of  the  angular  velocity  vector. 


p  +  qsm^tand-^r  cos^tan  0' 

e 

= 

qcos^-rsin^ 

(qsm<f’^-rcos^)  cosd 

Notice  that  for  small  angular  rotations,  as  expected. 


^  =  p: 

Q  =  q: 

V  =  r. 

At  this  point  the  kinematic  relationships  between  velocities,  as  seen  in  the  body  fixed 
frame,  and  the  rates  of  change  of  global  positions  and  Euler  angles  have  been  defined  The 
resulting  set  of  differential  equations  forms  a  consistent  set  in  that  given  a  set  of  vehicle 
velocity  data  versus  time,  its  position  and  attitude  may  be  computed 
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6.  Traoslational  Equations  Of  Motion 

The  global  acceleration  of  the  center  of  mass  is  derived  by  differentiating  the 
velocity  vector,  taking  into  account  that  the  center  of  mass  lies  in  a  rotating  reference 
frame  Considering  the  total  differentiation,  the  global  acceleration  of  the  center  of  mass 
becomes, 

^  =  v  +  djxpg +d)xd)xpj, +WXV,  (2  12) 

where  v  =  R„  The  translational  equation  of  motion  is  found  by  equating  this  acceleration 
to  the  net  sum  of  all  forces  acting  on  the  vehicle  in  three  degrees  of  freedom  (X,  Y,Z)  One 
important  factor  to  recognize  is  that  the  equation  of  motion  derived  in  this  manner  is  a 
vector  equation  with  the  components  expressed  in  the  body  fixed  frame  and  unit  vectors 
7.  j  and  k  This  has  been  deliberately  done  because  the  dominant  forces  acting  on  a 
submerged  body  in  motion  are  developed  in  relation  to  the  shape  of  the  vehicle  and  are 
more  conveniently  expressed  in  relation  to  the  body  axes  Equating  the  applied  force 
vector  to  the  acceleratioi.  results  in, 

r  =  m/^  +  tDxpj.+oix<bxpo+toxp£.^  (2  13) 

The  applied  force  vector  is  composed  of  gravitational  (weight)  and  hydrostatic  (buoyancy) 
forces  together  with  any  hydrodynamic  forces  arising  from  relative  motion  between  the 
vehicle  and  the  ocean  water  particles  These  will  be  discussed  in  further  detail  in  following 
sections 

7.  Rotational  Equations  Of  Motion 

To  develop  the  rotational  equations  of  motion,  the  sum  of  applied  moments 
about  the  vehicle's  center  of  mass  is  equated  to  the  rate  of  change  of  angular  momentum  of 
the  vehicle  about  its  center  of  mass  This  equating  will  provide  the  necessary  equations  of 
motion  for  the  remaining  three  degrees  of  freedom  In  the  practical  case  of  marine 
vehicles,  however,  the  statement  just  made  is  modified  slightly  because  it  is  much  more 


12 


difficult  to  assess  the  vehicle's  mass  moments  of  inertia  about  its  center  of  gravity  (CG) 
As  the  CG  changes  with  loading;  it  becomes  simpler  to  evaluate  the  mass  moments  of 
inertia  about  the  body  fixed  frame  which  lies  along  the  axes  of  symmetry  of  the  vehicle  in 
most  cases  The  mass  moments  of  inertia  for  the  vehicle  to  be  computed  are. 

'l  k  l' 

K=  L  4  4  (-14) 

.4  4  4. 

The  angular  momentum  of  the  body  is  thus, 

=  (2  15) 

resulting  in  the  total  applied  moments  about  the  origin  given  by, 

(2  16) 

Since  the  rate  of  change  of  angular  momentum  is  given  by, 

=  (2  17) 

and  the  acceleration  of  the  global  position  vector  is  given  by, 

4  =  ^  +  (2  18) 

then  the  rotational  equation  of  motion  in  vector  terms  is  given  by, 

=  /^w  +  djxC/^di^  +  m/pQ  X  v  +  Pq  xG)Xv}  (2  19) 

The  applied  moments  about  the  body  fixed  frame  arise  from  a  static  balance  of 
weight  and  buoyancy  effects  to  achieve  the  proper  trim  and  heel,  and  hydrodynamic 
moments  from  the  forces  applied  through  appendages  such  as  control  fins,  hydrodynamic 
effects  from  waves,  and  hydrodynamic  effects  from  relative  motion  between  the  vehicle 
and  the  water 

At  this  point,  there  are  three  translational  equations  obtained  from  equation 
2  13,  three  rotational  equations  obtained  from  equation  2  19,  and  six  unknown  velocities 
(diand  v)  In  itself  this  is  a  consistent  set  of  dynamic  equations  if  the  weight  and  buoyancy 


terms  are  always  self  canceling.  However,  with  weight  and  buoyancy  being  applied  in  a 
global  vertical  direction,  and  also  being  applied  at  different  locations,  they  represent  forces 
that  are  dependent  on  the  attitude  of  the  vehicle  Notice  that  without  weight  and  buoyancy 
forces,  the  six  consistent  equations  expressed  in  the  body  fixed  frame  of  the  vehicle  could 
be  solved  for  the  vehicle's  velocity  (given  applied  force  descriptions)  With  weight  and 
buoyancy  acting,  we  need  additional  equations  (constraints)  that  will  link  vehicle  attitude 
to  motion  so  that  the  combined  set  of  equations  can  be  solvable.  The  constraints  to  be  used 
are  developed  using  the  Euler  angles  which  are  a  general  set  of  angles  used  to  define  the 
attitude  of  the  rigid  body  From  these  angles  a  relationship  between  rate  of  change  of 
attitude  and  the  body  rotational  rates  given  earlier  as  the  angular  velocity  vector 
(a  =  f  p,q.r  J,  can  be  developed 

8.  Incorporation  Of  Vertical  Forces  Into  The  Equation  OfMotion 

The  weight  and  buoyant  forces  that  act  at  the  centers  of  gravity  and  buoyancy 
must  be  defined  from  static  analyses  For  submerged  bodies  the  weight  and  buoyancy  force 
vectors  do  not  change  with  vehicle  attitude  For  a  surface  ship,  B  will  change  with 
attitude,  and  the  righting  moments  must  be  computed  from  Naval  Architectural 
considerations  Assuming  that  weight  and  buoyancy  are  fixed  in  relation  to  the  body  fixed 
frame,  W  =  0l  +  0J  +  ( mg )K,  and  B  =  07  +  0J -(pV)K  Since  the  weight  and  buoyancy 
terms  in  the  applied  forces  act  in  the  global  vertical  direction,  they  must  be  transformed 
into  components  in  the  vehicle  fixed  frame  before  they  can  be  added  into  the  equations  of 
motion  Returning  to  equation  2  4,  we  therefore  find  the  components  along  the  vehicle 
fixed  frame  to  be  the  third  column  of  that  transformation  matrix  The  net  vertical  force 
components  become. 


-smB  ] 
7^  =(W-B)  cosQsm^ 
cosQcos^ 


(2  20) 


The  weight  portion  of  the  vertical  force  acts  at  the  center  of  gravity  of  the  vehicle  The 
buoyancy  portion  of  the  vertical  force  acts  at  the  center  ofbuoyancy  Because  these  forces 
act  at  positions  away  from  the  body  center  they  exhibit  a  moment  about  the  body  center 
given  by. 


-sm% 


-s/n0 


cosQsin^  ^  cos65/>t4> 

cos6cos<j>  CO50COS4> 


(221) 


This  moment  will  not  be  zero  even  if  W  and  B  are  identical  because  it  is  not  usually  the 
case  that  p^,  and  pg  are  collocated  For  static  stability  it  is  advisable  to  locate  the  center  of 
gravity  below  the  center  ofbuoyancy  The  total  vertical  force  vector  can  be  written  as. 


(2  22) 


and  is  added  negatively  to  the  left  hand  side  of  the  equations  of  motion 


9.  Development  Of  The  Full  Six  Degrees  Of  Freedom  Non-Linear  Equations 
Of  Motion  For  A  Marine  Vehicle 

Define  a  vector  x  of  vehicle  body  frame  velocities  to  be,  x  ■=  [u.v.w.p.qj  J , 
and  a  vector  z  of  global  positions  to  be  f  =  f  I ,  then  considering  M  as  a 

6x6  mass  matrix  including  translational  and  rotational  inertial  elements,  the  equations  of 


motion  can  be  written  in  the  following  vector  form, 

Mi  +  /(x)  +  F/z}  =  F^ 


(2  23) 


'■^g(x.z)  =  0 


(2  24) 


Therefore  with  suitable  knowledge  of  the  excitation  force  and  moment  loads, 
F^,  as  a  function  of  time  and/or  vessel  motion,  a  solution  for  the  vehicle's  dynamics  can  be 
obtained  A  more  detailed  insight  into  the  development  of  the  twelve  differential  equations. 


in  first  order  form  given  by  the  foregoing  analysis  shows, 

m^i"  +  ci)xp„/  +  m/a)X(oxpfj+djxv;  +  7g^2>=  Y,  .  (225) 

and. 

/  w  +  m/po  X  v/  +  djx//^to)  +  OT)^o  xtoxv/  +  m^f'r>=  Ki , 

It  helps  here  to  define  the  cross  product  coefficient  matrix  so  that, 
d)  X  Po  =  -Po  X  (0  =  -/ cros(p,5 ) /to 

where, 

‘o  ~yc, 

/cros(pG)y=  -z^  0  Xq  (2  28) 


(2  26) 


(2  27) 


Now  collecting  the  mass  matrix  coefficients  into  a  6x6  matrix  including  the  inertia  cross 
coupling  effects,  results  in  the  following  equation. 
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(2  29) 


The  remaining  terms  on  the  left  hand  side  of  the  equations  of  motion  arising  from  the 
centripetal  and  coriolis  accelerations  become. 


f(x)  = 


mfd)xcbxpc+0)xv^ 
dj X xGixv} 


(2.30) 


The  double  vector  cross  products  are  nonlinear  in  the  primary  velocity  variables  and  hence 
the  need  for  the  nonlinear  functional.  /(•)  The  reader  can  perform  the  indicated 
manipulations  to  express  individual  equations  within  the  set  if  so  desired 

The  components  of  the  hydrodynamic  and  external  forces  and  moments  acting 
on  the  vehicle  body  are  separated  in  the  above  analysis  into  six  components  each  acting 
along  the  vehicle  body  fixed  coordinate  axes  and  form  the  total  vector  of  forces  and 
moments  as, 

Fjt)  =  lX.(t),Yf(t).Z^(t).K,(t),M/t),N/l)]\  (2  31) 


where  the  vector  components  in  order  refer  to  the  surge,  sway,  heave  forces,  and  the  roll, 
pitch,  and  yaw  moments  respectively 

The  long  form  of  the  six  degrees  of  freedom  equations  of  motion  can  be  written 
as  follows. 


m[u-  vr  +  H(7  -  x^fq-  +r}  +  yj  pq-r)  +  zj  pr  +  qjJ  = 
-(iV -  BjsmQ+  X f 


(2  32) 


mfv  +  ur-wp  +  xjpq  +  r)- y^(  p'  +r  )+Za(qr  -  p)J  = 
{W  -  B)cosQsin^-^y^ 


(2  33) 


mfH-uq  +  vp  +  x^(pr-q)  +  y^(qr  +  p)-2jp‘+q^JJ=  - 

(2  34) 

(fV  -  B)cosQcos^  +  Zf 

hP  +  ( A  - +  L(P^ -^)-  JyJ q‘~r)-IJpq  +  r)-^ 

"»/ >G^'‘' -uq  +  vp)-2^fv  +  ur- wp)  J  =  ->'85;coseco5(j)  (2  35) 

-(2(;W -2QB)cosQsin^  +  Kj^ 

Kq  +  (l.-lJpr-IJqr  +  p)  +  IJ  pq-O+IJp'-r')- 

m[ x^fii'  -uq  +  vp)- ZqCu  -vr  +  wq ) ]  =  -( x^jW - XgB)cosQcos<^  (2  36) 

-{ ZgfV  -  2gB  JsinQ  +  S4f 

1/  +  (1^  - Ijpq - I^(p‘ -q‘)-Iy,(pr+q)+IJqr-p)  + 
m[xJv-\-ur~M/p)~y^(u-vr-irwp)]  =  (X(W-XgB)cos^sm^  (2  37) 

while  the  kinematic  relations  for  the  vehicle  rate  of  change  of  attitude  and  motion  over  the 
ocean  bottom  require  equations  2  11  and  2  8 

10.  Adaptation  Of  The  Non-Linear  Equations  Of  Motion  To  The  Vertical 
Plane 

Restricting  the  motions  of  the  vehicle  to  the  vertical  (dive)  plane,  the  only 
significant  motions  that  must  be  incorporated  to  effectively  model  the  vehicle  in  the  dive 
plane  are,  the  surge  velocity  (m),  the  heave  velocity  (w),  the  pitch  velocity  {q),  the  pitch 
angle  (0)  and  the  global  depth  positioii  (Z)  This  restriction  simplifies  the  twelve 
previously  developed  equations  to  a  system  of  four  non-linear  equations  of  motion,  which 
are; 

m(w-uq-2^q'  ~x^q)  =  Zf  (2  38) 
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- Tl - 

• 

I^q + mZf^wq  -  mx^  (w-uq)  =  Mf 

(2  39) 

• 

*) 

II 

(240) 

• 

Z  =  -usinQ  +  wcosQ 

(241) 

4r 

where. 

• 

Z,  =  Z^q  +  Z.vi  +  Z^uq  +  Zjtw  - 

Ip  f  C^(x)ip^^ix+(ff'-fi)cose+ii^Z,5,  +  Z^5j 

(242) 

• 

and, 

M^  =  M.q  +  M.w+Muq-¥M  uw 
f  q^  w  q  ^  w 

• 

1  xdx 

--pj  |k-X^| 

^  r«(/ 

(2  43) 

-( )cos  6  -  ( z^fV  -  z^B )  sin  0+  ( M^  S,  +  Mi  8^ ) 

• 

• 

results  from  expanding  the  Z^  and  the  Mf  terms  from  (2  31)  in  a  first  order  Taylor  series 

expansion  and  incorporating  both  hydrostatic  and  fluid  drag  forces 

• 

Equations  2  38,  2  39,  2.40  and  2  41  can  be  linearized  for  a  level  flight  path  when 

the  dive  plane  angle  is  zero,  i  e  5„  =  0.  By  setting  ail  time  derivatives  to  zero,  and 

neglecting  for  the  moment  the  hydrodynamic  drag  terms  the  following  are  obtained. 

• 

Zjiw  +  (W  -  B  )cosB  =  0 

(2  44) 

M^uw  -(XqW  -  XgB ) COS0  -(z^W  -  ZgBjsinQ  =  0 

(2  45) 

o 

II 

(2  46) 

• 

-«s/n0  +  uco50  =  0 

(2  47) 

• 

19 

• 

• 

•  ••••• 

• 

• 

• 

I 


If  the  assumption  is  made  that  the  vehicle  is  neutrally  buoyant,  then  it  can  be  said  that 


Xq  =  Xg  and  W  =  B  From  this  equations  2  44  and  2  45  can  be  reduced  to, 

=  0  (2.48) 

Mjuw-(Zfj-Zg)Bsm%-0,  (2.49) 

which  yield  the  nominal  position  =  0  and  si«0„  =  0,  resulting  in  the  solution  as 

either  =  0  or  =  Jt  Choosing  to  linearize  around  the  nominal  point  0^  =  0  results  in, 

=  =  0  (2  50) 

wq  =  w^q  +  q^w  =  0  (2.51) 

5/w0  =  cos0„0  =  0  (2  52) 

wcosQ  =  C-u'^sinQ^jB  +  fcosBgJw  =  w  (2.53) 

The  linear  equations  of  motion  are  then  written  as; 

mXf.  )q  =  Zjm  -^(Z^  +  m  )uq  +  ZgU‘h  (2  54) 

( I ^  -  M^)q-(M^+  mxg  )w  =  A/,  ww  +  (M^-  mx^  )uq  -(zQ-Zg)W^  +  Mgurb  (2  55) 

0  =  <7  (2  56) 

Z  =  -ue+w  (2  57) 


As  equations  2  42  and  2  43  show,  both  Zg  and  A/§  are  a  linear  combination  of 
the  respective  stem  and  bow  hydrodynamic  control  surface  coefficients  and  the  respective 
input  value  of  5  This  makes  the  system  of  equations  as  a  multiple  input  system  To 
reduce  this  system  into  a  single  input  system  the  linear  combination  of  control  inputs  will 
be  modified  into  the  following  form, 

Z*=rZj  +cxZj  )  (2  58) 

s  s 

This  will  allow  a  single  input  6  to  control  both  stem  planes  and  bow  planes,  and  will 
cause  the  bow  planes  to  be  slaved  to  the  stem  planes  This  technique  is  known  as  dual 
control  The  value  OC  will  range  from  - 1  to  1  The  selection  of  the  value  of  Ot  will  allow 
the  planes  to  operate  as  desired  for  the  particular  maneuvering  condition,  i  e.,  a  =  0  for  no 
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bow  plane  control,  a  =  -1  for  bow  plane  and  stem  plane  control  opposed  to  each  other, 
yielding  the  maximum  pitch  moment,  and  a  =  1  for  bow  and  stem  plane  control  in  the 
same  direction,  yielding  the  maximum  heave  force 

In  state  space  form  these  equations  can  be  written  as, 

'el  r  0  0  1  oTel  r  0  ■ 

w  ^  a,,"  ®  ^  5  (2  59) 

q  cum  0  q  b-M' 


[z\  L  -«/  10  olzj  L  0  J 

where  the  coefficients  and  6,  are  given  by; 

+  (2  60) 

a„A  =  (ly  -  M^)Z^  +  (mx^  (2  61) 

a,,D,  =  (I^  -  M^Xfn->r  Z^)  +  (mxa  +  Z^X -mx^)  (2  62) 

a^^D^  =  -(niXQ->rZ^)W  (2  63) 

bA  =  (ly  -  M^)Zi  +  (mxo  +  Z^)M,  (2.64) 

a, ^  =  (m-Z^)M^  +  (mX(;  +  M^)Z^  (2  65) 

a,,D,  =(m-Z^XM<,-mXi;)+(mxQ  +  M^Xfff^Z^)  (2  66) 

a,,D^  =  -(m-ZJW  (2  67) 

b, D^=(m-Z^)M^  +  (mx^  +  M,)Z^  (2  68) 

and  z^g  =  Z(y-Zg  is  the  metacentric  height 


This  linear  system  of  equations  becomes  the  basis  for  the  control  law  design  in 
the  following  section 

B.  CONTROL  LAW  DESIGN 
1.  Introduction 


The  control  design  problem  can  be  stated  as  follows:  Given  the  system. 


x  =  Ax-^Bh, 


(269) 


where  the  state  vector  equation  is. 


6 

w 


X  = 


<7 

Z 


(2  70) 


how  do  we  find  5,  such  that  the  system  will  behave  as  desired.  The  type  of  control  that  is 
of  interest  in  this  problem  is  closed  loop  control,  where  5  is  a  function  of  the  state  x 
Since  the  state  x  is  used  to  determine  the  control  effort  5('x^  it  is  called  feedback  control. 

2.  Pole  Placement 

A  linear  full  state  feedback  control  law  is  introduced  in  the  form 


b  =  -Kx,  (271) 

where  K  is  the  feedback  gain  vector  to  be  determined  such  that  the  closed  loop  system  of 
equations  2  69  and  2  71  has  the  desired  system  dynantics.  Substituting  equation  2  71  into 
equation  2  69  yields, 

x  =  (A-BK)x  (2  72) 


The  actual  characteristic  equation  of  this  closed  loop  system  is  given  by 

del[A-BK-sl]^0  (2  73) 

The  gain  vector  K  can  be  chosen  such  that  the  actual  characteristic  equation  assumes  any 
desired  set  of  eigenvalues  If  the  desired  locations  of  the  closed  loop  poles  are  chosen  at 
5  =  5,  for  /=1,  ,n,  the  desired  characteristic  equation  becomes 

<5-5|/5-5-/.Y5-5,>  =  0  (2  74) 

The  required  values  of  K  are  obtained  by  matching  coefficients  in  the  two  polynomials  of 
the  actual  and  desired  characteristic  equations 
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Consider  the  previously  developed  linear  state  matrix  equation  2  59 


♦ 


Substituting  equation  2  71  into  equation  2  59  results  in  the  closed  loop  system 
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The  characteristic  equation  of  the  closed  loop  system  is  given  by 


(2  75) 


j  -5  0  10 

*^t3“Gs  a^iU-b^u~k,-s  a^M-httk,  -hu'k^ 

det\  ,  ‘  , 

Ks'ca  ~  ^“‘^1  ~  b-M^k,  a^u  -  b-ju'k.^  -  s  -b,u'k^ 

\  -U  1  0-5 

which  after  some  algebra  reduces  to 

5'*  +  ^  A.k,  +  A^k^  -  £,  )s^  ■¥(-B^k^-  B^k.,  -  B-^k^  -  B^k^  -  £,  ).r 
+  (-C,k,  -  C;k;  -  C^k^  -  Eyjs+(-D^k^  -  D,k^)  =  0 

where , 

A  =  -B^  =  b^w 

A,  =  -B^  =  b,u' 

B,  =(a.J}^  -a^,b,  Ju^ 

C,  —  Dy  —  ^^i-<ib\  ~  ib,  )ZQgU 

C,  =  +£  -a,,E)u^ 

a  =  ra,A-a,,A,V 

£,  =  ^'a,,  +a.,Ju 


(2  76) 


(2  77) 


(2  78) 
(2  79) 
(2  80) 
(2  81) 
(2  82) 
(2  83) 
(2  84) 
(2  85) 
(2  86) 
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^w^zjJ^gb^  (2  87) 

Now  if  the  closed  loop  poles  are  placed  at  -p,,  -p,,  -p,.  and  -p^,  then  the 
desired  characteristic  equation  is, 

{s  +  pj(s  +  p,)(s  +  p^Xs  +  pj  =  0,  (2  88) 


or 


5^  +  0,5^  +  aS'  +  a,5+ =  0 


(2  89) 


with, 

«i  =Pi  +  P:-*-A  +  A  (2  90) 

«:  =P^P.+  PPi  +  AA  +  A  A  +  PiP*  +  A  A  (291) 

=  PPiPi  +  PPzP,  +  PiPJ>*  +  A: AP4  (2  92) 

«4=AA:AA4  (2  93) 

The  control  gains  can  now  be  computed  by  equating  coefficients  of  the  actual  and  desired 
characteristic  equations. 


Ak,  +  =  -a,  -  £,  (2  94) 

+  B^k,  +  +  B^k^  =  a,  +  £;  (2  95) 

r,^,  +  C£+CA=a3  +  ^,  (296) 

(D\  +  D2Jki=a4  (2  97) 


Now  that  a  method  for  computing  the  gains  of  a  controllable  single  input  system 
that  will  place  the  poles  at  any  desired  location  has  been  developed,  the  selection  of  pole 
location  must  be  accomplished 

3.  Pole  Location  Selection 

In  a  typical  second  order  system  control  law  design  the  transient  response 
specifications  will  be  given  This  results  in  an  allowable  region  in  the  s-plane  where  the 
desired  location  of  the  poles  can  be  obtained  For  higher  order  systems  the  concept  of 
dominant  roots  can  be  employed  In  selecting  poles  for  a  physical  system  it  is  necessary  to 
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look  at  the  physics  of  the  system  If  the  poles  are  specified  too  negative,  a  very  small  time 
constant  for  the  control  system  will  result,  and  the  physical  system  may  not  be  able  to 
react  that  fast  The  control  law  u  =  -Kx ,  implies  that  for  a  given  state  x  the  larger  the 
gain,  the  larger  the  required  control  input  In  practice  there  are  limits  typically  placed  on  u 
(actuator  size  and  saturation)  Occasional  control  saturation  is  not  serious  and  may  even 
be  desired  A  system  that  never  saturates  is  in  all  likelihood  over  designed 

Considering  the  control  law  design  to  stabilize  the  submarine  to  a  level  flight 
path  at  6  =  0  it  will  be  required  that  the  submarine  return  to  level  flight,  after  a  small 
disturbance  in  0  or  Z,  within  the  time  it  takes  for  the  vehicle  to  travel  three  ship  lengths 
Since  the  submarine  is  about  14  feet  long  and  it  travels  at  5  ft/sec,  the  required  recovery 
time  is  about  10  seconds  This  means  that  the  time  constant  is  about  3  seconds,  and  the 
closed  loop  poles  should  be  placed  at  approximately  -0. 3 

Control  design  using  pole  placement  is  very  easy  using  MATLAB  The 
appropriate  MATLAB  command  is  place,  which  accepts  as  inputs  the  A  and  B  matrices 
along  with  a  vector  of  the  desired  closed  loop  poles,  and  returns  the  gain  vector  K.  Using 
equation  2  59  with  a  nominal  speed  of  5  ft/sec  and  the  vector 
p  =  [-0.3  -0.31  -0.32  -0.33/,  (place  does  not  like  poles  in  the  exact  same  location), 
the  gain  vector  K  was  calculated  using  MA  TLAB  and  by  simultaneously  solving  equations 
2  94  through  2  97  yielding  the  same  results  Substituting  the  gain  vector  and  state  variable 
vector  into  equation  2  71  results  in  a  control  law  of 

5  =  -r-O.99170-O.8333K-O.6O26^  +  O.O351Z/  (2  98) 

Using  this  control  law  along  with  equations  2  38  through  2  43  a  simulation 
program  was  developed  to  investigate  the  vehicle's  response  to  initial  disturbances,  the 
results.of  which  will  be  presented  in  the  following  chapter 
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in.  PROBLEM  IDENTIFICATION 


A.  VEHICLE  SIMULATIONS 

With  the  control  law  developed  and  the  vertical  plane  equations  of  motion  defined,  a 
vehicle  simulation  program  was  developed  (see  Appendix  A)  This  program  was  used  to 
simulate  the  vehicle's  response  at  a  variety  of  speeds  and  initial  disturbances  The 
simulations  were  conducted  with  several  simplifications  applied  to  the  vertical  plane 
equations  of  motion  (E  O  M ),  equations  2  38  through  2  43  The  simplifications  or 
constraints  applied  to  the  E  O  M  were  , 

1 )  to  consider  the  body  drag  forces  negligible, 

2)  to  consider  the  vehicle  to  be  neutrally  buoyant,  i  e  ,  ( W=B), 

3)  to  assume  that  the  locations  of  the  center  of  gravity  and  center  of  buoyancy, 
in  the  X-Y  plane,  are  coincidental,  and 

4)  that  the  rudder  is  restricted  to  ±  23  degrees 

These  assumption  resulted  in  the  reduction  of  equations  2  38  through  2  43  to  the 
following  system  of  equations; 

10  ®  r  ^ 

0  (m-Z,)  -Z^  0  \v  m(uq  +  Z(.q^ J  +  Z^uq  +  Z^uw  +  u'‘Zi 

0  -M.  I^-M^  0  q  m(-ZQwq)  +  M^uq+  M^uw-  Bz^gSinQ  +  u'Mi^ 

0  0  0  1  Z  -usind  +  wcosQ 

forming  the  basis  of  the  initial  vehicle  simulations 


The  vehicle's  response  to  small  disturbances,  as  stated  earlier,  was  investigated  over  a 
wide  range  of  speeds  A  small  sample  of  these  simulations  is  shown  in  Figures  3-1  through 
3-3 

As  Figure  3-1  depicts,  the  vehicle  returns  to  level  flight  in  approximately  10  seconds. 
This  response  is  as  expected  based  on  the  control  system  design  of  Chapter  II  The 
simulation  results  represented  in  Figure  3-2  also  show  the  vehicle  steadying  out  in  a  level 
flight  at  the  desired  depth,  however  the  vehicle  requires  a  significant  amount  of  time  for 
this  to  occur  The  reason  for  this  extended  amount  of  time  is  that  the  control  law  gains 
were  set  for  the  nominal  operating  speed  of  5  fps,  and  the  vehicle  is  being  simulated  at 
one-third  of  that  speed  This  results  in  a  very  slow  return  to  ordered  position  caused  by 
the  reduced  hydrodynamic  effects  at  this  speed 

The  final  plot  of  the  sample  vehicle  simulations.  Figure  3-3,  illustrates  a  problem  with 
the  vehicle  stability  As  can  be  seen,  the  vehicle  does  not  return  to  the  desired  position,  but 
instead  steadies  out  at  a  pitch  angle  of  eight  degrees,  with  a  steady  state  dive  plane  value 
of  23  degrees  This  same  type  of  response,  but  with  different  pitch  angle  values,  was 
observed  in  all  the  simulations  conducted  below  approximately  19  fps  These  results 
indicate  that  there  is  a  critical  point,  or  speed  at  which  the  vehicle  stability  is  affected. 
Determining  this  point  will  be  the  basis  for  the  remainder  of  this  chapter 

B.  CRITICAL  SPEED  IDENTIFICATION 
1.  Eigenvalue  Analysis 

Consider  the  nonlinear  system  of  state  equations  3  1,  which  can  be  written  in  the 

form, 

x  =  f(x)  (3  2) 

It  is  known  that  the  equilibrium  points,  of  the  system  are  defined  by 
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non  dimensional  depih  theta/dcita  (degrees) 


15 


THETA/DELTA  vs.  TIME 


1.2 


non-dimeiuional  time 
DEPTH  vs.  TIME 


0.8- 


0.6  - 


Figure  3-1.  Vehicle  Response  For  a  Nominal  Operating  Speed  of  5  fps. 
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/rxj  =  0 


(3.3) 


This  is  a  nonlinear  system  of  algebraic  equations  and  it  may  have  multiple  solutions  in 
which  means  that  the  nonlinear  system  of  equations  may  have  more  than  one  position  of 
static  equilibrium  If  one  equilibrium  value  x^.  is  selected,  the  stability  properties  can  be 
established  by  linearization  The  linearization  converts  equation  3  2  into  the  form 

x=Ax.  (3  4) 

where  A  is  the  Jacobian  matrix  of /(x)  evaluated  at  the  equilibrium  point  x^. 


A  =  ^M 

Bxl 


(3  5) 


and  the  state  x  has  been  redefined  to  designate  small  deviations  from  the  equilibrium 

x-4x-x„  (36) 

As  long  as  all  eigenvalues  of  A  have  negative  real  parts,  the  linear  system  will  be  stable 
This  means  that  the  equilibrium  x^  will  be  stable  for  the  nonlinear  system  as  well  This 
statement  is  nothing  more  than  Lyapunov's  linearization  technique 

The  question  that  must  be  answered  is,  what  happens  if  one  real  eigenvalue  of 
the  linearized  A  matrix  is  zero'’  The  interesting  case  here  is  when  the  rest  of  the 
eigenvalues  have  all  negative  real  parts,  otherwise  x„  is  unstable  and  the  problem  is  solved 
If  the  case  of  a  zero  eigenvalue  appears  too  specialized  to  be  of  any  practical  use,  consider 
this  Assume  that  /(x)  depends  on  one  physical  parameter  and  that  physical  parameter  is 
allowed  to  vary  over  some  range  Then  clearly  A  will  depend  on  that  parameter  and  as  the 
parameter  varies,  it  is  possible  that  one  real  eigenvalue  of  A  will  become  zero  for  a  specific 
value  of  the  parameter  The  problem  is  then  to  establish  the  dynamics  of  the  nonlinear 
system  as  one  real  eigenvalue  of  A  crosses  zero,  i  e  ,  goes  from  negative  to  positive  As 
the  solution  evolves  in  time,  things  are  interesting  only  along  the  direction  of  the 
eigenvector  that  corresponds  to  the  critical  eigenvalue  Along  the  rest  of  the  directions  in 
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the  state  space,  everything  should  converge  back  to  the  equilibrium,  remember  it  was 
assumed  that  all  remaining  eigenvalues  of  A  have  real  negative  parts  Although,  strictly 
speaking,  this  is  a  true  statement  for  linear  systems,  there  are  technical  reasons  that  force  it 
to  be  true  for  nonlinear  systems  as  well,  the  only  difference  is  that  the  corresponding 
directions  in  the  state  space  are  curved  instead  of  straight. 

By  taking  the  previously  developed  characteristic  equation  2  77,  and 
remembering  that  the  roots  of  this  equation  are  the  eigenvalues  of  interest,  it  is  easy  to  see 
that  if  of  equation  2.93,  is  set  equal  to  zero,  then  one  eigenvalue  will  be  zero  If  this  is 
done,  and  recalling  that  holds  a  non-zero  value  equation  2  97  reduces  to 

^^1 1^2  ^gb( ~  ^ ■  (3  7) 

Since  all  parameters  of  equation  3  7  have  a  fixed  value  with  the  exception  of  «,  there  must 
be  some  value  of  u  that  causes  the  linear  A  matrix  to  become  unstable  Recognizing  this 
fact  will  allow  equation  3  7  to  be  solved  in  terms  of  u  yielding 


y  -  ^08 ^^2 A 


Using  equations  2  60  through  2  68,  equation  3  8  can  be  simplified  into  the  following  form, 

r  ni  : 

u=  - HoihE -  (3  9) 

Substituting  the  appropriate  values  from  Appendix  B,  equation  3  9  can  be  solved  for  the 
value  of  M  that  causes  the  systems  to  become  unstable  This  results  in  a  value  of  w  = 

1  8979  fps  using  an  a  =  0  This  value  corresponds  very  closely  to  the  value  of  1  9  fps 
obtained  by  vehicle  simulations 
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As  a  check  to  this  solution,  a  MATLAB  program  was  written  to  compute  the 
closed  loop  eigenvalues  of  the  matrix  equation  2.59  as  a  function  of  the  speed  u 
Figure  3-4  shows  a  plot  of  the  real  parts  of  the  eigenvalues  versus  speed  It  can  be  seen 
that  the  real  part  of  the  eigenvalue,  represented  by  the  dash-dot  line,  becomes  positive  at 
speeds  less  than  approximately  19  fps  supporting  the  earlier  findings  (For  clarity 
purposes,  the  values  of  the  eigenvalue  represented  by  the  dash-dot  line  were  scaled  by  a 
factor  of  ten  ) 

2.  Steady  State  Analysis 

To  determine  the  existence  of  multiple  steady  state  solutions,  for  the  system 
represented  by  equation  2  59,  it  is  necessary  to  perform  a  steady  state  analysis  on  the 
system  of  equations  which  models  the  vehicle  dynamics  Using  matrix  equation  3  1.  and 
setting  all  time  derivatives  equal  to  zero,  results  in  the  following  set  of  equations. 


q  =  0,  (3  10) 

Z^uq  +  Z^uw  +  tt-Zj6  =  mf-uq-Zaq- (311) 
M^uq  +  A/.ww  -  Bz^g  sinB  +  u'Mgb  =  ntz^wq ,  (3  12) 

-«  5/>f0  +  wcosG  =  0  (3  13) 


Simplifying  these  four  equations  yields  the  following  two  equations  in  terms  of  6 ,  u.  and 
the  physical  parameters  of  the  vehicle, 

Z^u' tanQ  +  u'Zg5  =  0  (3  14) 

M^u' lanQ~BZ(jgSinQ  +  u‘Mgb  =  0  (3  15) 

Now  if  equation  3  14  is  multiplied  by  Mg,  equation  3  15  is  multiplied  by  Zj  and  the  two 
resulting  equations  are  set  equal  to  each  other  and  simplified,  the  following  single  equation 
is  obtained. 

/ ^ Z,Mj  -  AY.Zj yM‘  +  ZgBZfjg cosQ / sinQ  =  0  (316) 
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Through  inspection  it  can  be  seen  that  equation  3.16  could  have  multiple 
solutions  in  0 ,  besides  the  trivial  solution  of  0  =0  So,  if  this  equation  is  rearranged  and 
solved  for  0*0,  solutions  occur  when 


CO50  = 


(3  17) 


Therefore,  the  steady  state  value  of  0  is  represented  by 


0.,  =  cos 


(3  18) 


which  may  have  multiple  solutions  based  on  the  value  of  the  speed  u  As  discussed 
previously,  these  multiple  solutions  were  evident  in  the  simulations  conducted  earlier  in 
this  chapter 

Knowing  that  the  maximum  value  for  co50  is  equal  to  one,  the  right-hand  side 
of  equation  3  17  will  be  true  for  values  less  than  or  equal  to  one  If  this  constraint  is 
imposed  on  3  17,  then  the  equation  may  be  manipulated  into  the  following  form 


«-  < 


^8^*08 


(M,Z,-Z,M,) 


(3  19) 


Rearranging  this  equation  and  solving  for  the  upper  limiting  case  yields  the  same  results  as 
equation  3  9  Again  this  supports  the  critical  speed  values  discussed  earlier 

Equation  3  19  can  be  expressed  in  a  non-dimensional  Froude  number  form  by 
dividing  the  equation  by  both  the  gravitational  constant  g  and  the  physical  parameter  zgg 
and  taking  the  square  root  of  the  result  This  will  convert  equation  3  19  into  the  non- 
dimensional  form  given  below. 


Fn  = 


GB 


Z,B 


g(M^Z^-Z,Mi,) 


(3  20) 
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With  the  steady  state  6  solution  derived,  equations  for  both  the  steady  state 
value  of  5  and  steady  state  value  of  Z  can  be  determined.  Using  equation  3  14  and 
substituting  the  steady  state  value  of  6  will  allow  the  steady  state  5  equation  to  be 
obtained, 

(3  21) 

A 

To  achieve  the  steady  state  solution  of  Z.  begin  by  writing  equation  2  71  in 
expanded  general  form  as 

8  = -{ k^Q  +  +  k^q  +  k^Z  J  (3  22) 

Then  by  applying  the  same  conditions  to  equation  3  22  that  resulted  in  equations  3  14  and 
3  15,  and  using  the  previously  developed  steady  state  equations  for  6  and  5,  an  equation 
of  the  form 

+  k.u  tanQ^ 

Z„=— 5 - ^ -  (3.23) 

-K 

is  obtained 

The  loss  of  stability  of  an  equilibrium  and  the  generation  of  additional 
equilibrium  states,  is  called  a  pitchfork  bifurcation  and  is  very  common  in  nature  The 
buckling  of  a  beam  is  one  such  example  Using  equations  3.18,  3  21,  and  3  23  as  a  basis 
for  another  MATLAB  program,  the  steady  state  solutions  of  0,  5,  and  Z  versus  Froude 
number  were  investigated  with  the  results  displayed  in  Figures  3-5  through  3-7  These 
three  plots  are  referred  to  as  supercritical  pitchforks,  so  named  because  upon  the  loss  of 
stability  of  the  trivial  equilibrium  the  additional  nearby  equilibrium  states  are  stable 
Graphically,  this  can  be  represented  in  Figures  3-5  through  3-7,  where  the  solid  curves 
represent  stable  and  dotted  curves  represent  unstable  equilibria  [Ref  17]  These  three 
figures  also  demonstrate  the  control  input  saturation  Upon  control  surface  saturation,  the 
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steady  state  equations  for  9 ,  5  and  Z  are  no  longer  valid  due  to  the  limit  placed  on  the 
maximum  plane  angle  The  steady  state  equations  that  hold  for  this  region  of  saturation 
are  given  below. 


8„  =  ±  22.9°  (0  4  radians ). 


2 


(3  24) 
(3  25) 


and. 


0„  =  sin 


(3  26) 


As  can  be  seen  in  Figures  3-5  and  3-6,  at  an  approximate  Froude  number  of  1  05  the 
control  surfaces  saturate  resulting  in  the  linear  solutions  displayed  below  that  saturation 
Froude  number  In  Figure  3-7  there  is  no  steady  state  solution  because  if  the  steady  state 
values  of  6  and  6  are  placed  into  the  Z  equation,  Z  is  non-zero  below  the  saturation 
Froude  number 

A  parameter  introduced  in  Chapter  II  was  a,  which  allowed  us  to  go  from  a 
multiple  input  system  to  a  single  input  system  This  parameter  has  a  significant  effect  on 
the  location  of  the  critical  speed  and/or  Froude  number  Figure  3-8  shows  the  relationship 
been  the  critical  speed  and  the  metacentric  height  zgg  for  different  values  of  a 

These  three  curves  shown  in  Figure  3-8  can  be  converted  into  a  single  curve  by 
plotting  a  versus  Froude  number,  and  is  shown  in  Figure  3-9.  As  can  be  seen  in  this  plot 
the  critical  Froude  number  varies  from  0  81  to  1  22  over  the  allowable  range  of  a 

In  addition  to  the  effect  that  a  has  on  the  location  of  the  critical  Froude  number, 
it  also  effects  the  magnitude  of  the  steady  state  values  at  a  given  Froude  number  and  the 
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STEADY  STATE  VALUES  OF  THETA  vs.  FROUDE  NUMBER 
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Figure  3-5.  Steady  State  Values  of  6  vs  Froude  Number. 
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STEADY  STATE  VALUES  OF  DELTA  vs.  FROUDE  NUMBER 
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Figure  3-6.  Steady  State  Values  of  5  vs  Froude  Number. 


39 


Z  (feel) 


(sclj)  jnn 


15 


•»  _ 


1.5- 


0.5 


( 


CRITICAL  SPEED  vs.  METACENTRIC  HEIGHT 
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Figure  3-8.  Relationship  Between  Critical  Speed  and  Metacentric  Height  For 

Various  Values  of  a. 
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Figure  3-9.  Froude  Number  as  a  Function  of  a. 
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location  at  which  the  control  surfaces  saturate  When  comparisons  are  made,  it  can  be 
shown  that  as  the  value  of  a  becomes  more  positive  the  location  of  the  critical  Froude 
number  and  steady  state  values  of  6  and  Z  become  greater  Comparing  the  changes  in  the 
5  steady  state  solutions,  the  shape  of  the  5  pitchfork  does  not  change  however,  the 
location  of  the  critical  point  moves  in  the  direction  of  increasing  Froude  number  as  a 
increases  The  results  of  this  discussion  are  shown  graphically  in  Figures  3-10  through 
3-12 

An  interesting  result  is  displayed  in  Figure  3-13  This  figure  depicts  the 
relationship  between  a  and  the  critical  or  saturation  Froude  number  At  the  lower  Froude 
numbers,  there  is  very  little  difference  between  the  critical  Froude  number  and  the 
saturation  Froude  number  This  demonstrates  that  upon  reaching  a  low  critical  Froude 
number  the  control  surface  immediately  saturates  attempting  to  keep  the  vehicle  stable  At 
the  higher  Froude  numbers  there  is  a  significant  difference  between  the  critical  Froude 
number  and  the  saturation  Froude  number  This  is  because  as  the  Froude  number 


approaches  2  78,  the  value  used  in  the  control  law  design  conducted  in  Chapter  II,  the 
control  surfaces  will  not  have  to  saturate  to  maintain  stability 

3.  Controllability  Analysis 


In  general,  every  system  of  the  form. 


x=  Ax-^  Bu. 
>=Cx, 


(3  27) 


can  be  divided  through  a  series  of  transformations  into  four  subsystems 

1 )  A  controllable  and  observable  part 

2)  An  uncontrollable  and  observable  part 

3 )  A  controllable  and  unobservable  part 

4)  An  uncontrollable  and  unobservable  part 
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STEADY  STATE  VALUES  OF  THETA  vs.  FROUDE  NUMBER 


-  ibr-d^pha  =  0 


-  for  alpha  =  1 


for  metacentric  height  0.1  ft 


Figure  3-10.  Comparison  of  Steady  State  Values  of  6  for  Different  Values  of  a. 
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STEADY  STATE  VALLTS  OF  DELTA  vs.  FROLDE  NUMBER 
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Figure  3-11.  Comparison  of  Steady  State  Values  of  5  for  Different  Values  of  a. 
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Figure  3-12.  Comparison  of  Steady  State  Values  of  Z  for  Different  Values  of  a. 
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Figure  3-13.  Relationship  Between  Critical/Saturation  Froude  Number  and  a. 
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This  is  known  as  Kalman's  decomposition  theorem.  Now,  since  the  transfer  hmction  of 
any  system  is  determined  by  the  controllable  and  observable  subsystems  only,  then  the 
transfer  function  may  contain  less  information  than  is  actually  needed  to  model  the 
complete  system  The  precise  definition  of  controllability  is: 

A  system  is  said  to  be  controllable  if  any  initial  state  x(t^)  can  be  driven  to  any 
final  state  x(tj)  using  possibly  unbounded  control  u(t)  in  finite  time  /„<  t  •  tf 
From  the  state  equations  3  27,  this  should  depend  only  on  A  and  B,  where  A  is  i  n  x  n 
square  matrix  The  criterion  for  controllability  is  as  follows:  Compute  the  controllability 
matrix 

c  =  [B.AB,A^B . A'^'B],  (3.28) 

and  the  system  is  controllable  if  and  only  if  the  rank  of  c  (the  number  of  linearly 
independent  rows  or  columns)  is  n.  Roughly  speaking,  c  shows  how  possible  it  is  to 
change  the  state  of  a  system  using  the  input.  For  a  single  input  system  B  xsn  x  1  and  c  is  a 
square  matrix  The  test  is  then  that  c  be  nonsingular 

det(c)*0.  (3  29) 

Now.  if  the  controllability  matrix  is  formed  using  the  A  and  B  matrices  from 
equation  2  59  then  the  following  matrix  results. 


where. 


^13 

^^14 

c,. 

t^:4 

<^33 

^34 

^4: 

<^43 

=  0. 

=  A.u" 

=  0,,M 

a,Jb 

(3  30) 


(331) 
(3  32) 
(3  33) 


•  • 


c,j  =by-(a^^a,y  +a,^a.,u'J  +  b,u~(a.jZ(jg+a^^a^,u'  +aJu-J, 

(3  34) 

•J 

c,,  =by\ 

(3  35) 

• 

c,;  =a^^u^b^ 

(3  36) 

Ar 

c,,  =by-(a^^'u-  +a,^u'  )-^b,u‘(a^j2Qg  +a,,ai,u-  +a,M„u'), 

(3  37) 

• 

6;  w  Va,  j^oga, , « + a,  + o^.u(a^ ,  ’w'  +  a,  ,a,|  u*  >  -i- 

(3  38) 

c„  =ku\ 

(3  39) 

• 

Cj,  =ayU^bf  +a.,u^b,. 

(3  40) 

c„  =b,u'{'a,,a,,u' +a,ja,,u'J  +  b.u'^a,jZ^g+a,ia,.u'  +a,,'w'>. 

(341) 

c„  =6,MVa,|M("<3||a,|«‘  +a|,a,,M‘  +a,,‘w-»  + 

b.u'(a^  jZ(jga,^u  +  a,,a,yZQgU  +  a^^u( a^ )  + 
a.Ma,^ZQg  +a|,a,,«*  +a„~u' )), 

(3  42) 

• 

o 

II 

(3  43) 

c,.,  =by-. 

(3  44) 

• 

• 

=a^yb^  +(a^,u-u)Lu‘, 

(3  45) 

and  =b^u'('a^^'u'  +a2^u(a^,u-u))  +  b,u‘(a^yZ(;g^^■a^^a^,u^  +a,.u(a^,u-uj) 

(3.46) 

Now,  if  the  determinant  of  the  controllability  matrix  is  computed  as  a  function  of 
Froude  number,  then  it  can  be  shown  that  the  system  becomes  uncontrollable  at  the  critical 
Froude  number  The  graphical  results  of  this  computation  are  shown  in  Figure  3-14  for 
different  values  of  the  parameter  a 

With  the  critical  parameter  that  causes  the  vehicle  to  lose  stability  identified,  i  e  , 
the  Froude  number,  and  the  effect  that  the  parameter  a  has  on  the  location  of  this  critical 
value  understood,  the  simplifications  discussed  at  the  beginning  of  this  chapter  can  be 
removed  from  the  E  O  M  and  an  analysis  on  the  effects  of  incorporating  these  additional 
parameters  can  be  conducted  This  detailed  analysis  will  be  the  topic  of  Chapter  IV 
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Figure  3-14.  System  Controllability  vs  Froude  Number. 


50 


IV.  BIFURCATION  ANALYSIS 


A.  ASYMMETRIC  PITCHFORK 

Now  that  the  critical  parameter,  the  Froude  number,  has  been  identified,  the 
simplifications  applied  to  the  E  O  M  in  the  previous  chapter  can  be  removed  The 
resulting  system  of  equations  is  given  by  equations  2  38  through  2  43  To  determine  the 
existence  of  multiple  steady  state  solutions  for  this  new  system  of  equations  it  is.  once 
again,  necessary  to  perform  a  steady  state  analysis  as  conducted  in  Section  B  2  of  the 
previous  chapter  This  analysis  will  reduce  the  system  of  equations  into  the  following, 

Zjftv  -  ^  pQ/4h'|h'|  +  JcosB  +  ZjM'S  =  0  (41) 

-XgBjcos6-/'Z(jfy-:gBMfnd-h  M^u'5  =  0,  (4  2) 

with  equation  3  13  also  obtained  from  the  steady  state  Z  equation  By  multiplying  equation 
4  1  by  A/j  and  equation  4  2  by  Zj  and  setting  the  resulting  equations  equal  to  each  other 

the  following  equation  is  obtained; 


(4  3) 


CCtV-BjMg+('X(jtV-XgBjZ^J  cosQ  +  5/rt  6  =  0 

This  equation,  along  with. 

w'  =  «/a«0,  (4  4) 

from  equation  3  13  is  the  exact  steady  state  solution  set  with  all  parameters  included 


Figure  4-1  shows  the  results  obtained  when  equation  4  3  is  solved  numerically,  using  the 
FORTRAN  program  in  Appendix  A  for  a  small  value  of  (the  difference  in  the  location 
of  the  center  of  gravity  and  buoyancy),  with  the  submarine  neutrally  buoyant  = 
This  small  value  of  perturbs  the  pitchfork,  as  seen,  from  the  symmetric  pitchfork 
displayed  earlier  in  Figure  3-5  This  comparison  shows  that  the  critical  or  bifurcation  point 
has  moved  to  a  lower  Froude  number,  with  the  stable  solutions  represented  by  the  solid 
lines  and  the  unstable  solutions  represented  by  the  dashed  lines  in  both  cases.  The  stable 
steady  state  solution,  that  the  vehicle  will  acquire,  is  dependent  on  the  initial  conditions 
that  are  given  to  the  vehicle  Using  the  results  displayed  in  Figure  4-1  as  an  indication  that 
may  also  be  a  critical  parameter,  the  subsequent  task  is  to  determine  the  relationship 
between  the  critical  Froude  number  and  xq^  or  /f  Fn.x^g) 


B.  BIFURCATION  GRAPHS 

The  first  order  of  business  in  the  attempt  to  identify  f(Fn.Xgg)  is  to  simplify  the 
exact  pitchfork  equation  This  can  be  accomplished  by  replacing  the  co59  and  5W0  terms 
in  equation  4  3  with  the  following  Taylor  series  expansions 

w  I  w’  2wu'-w- 

sinQ  = - - = - — ^ -  (4  5) 

u  2  u  2u 

cosQ  =  \-  —  ^  =  — — (4  6) 
2  2m- 

which  will  allow  the  exact  solution  to  be  represented  as  follows. 


Z(,(ZaiV-:gB)w^  +  2u^CcA{  Mg-x^ZgJw\w\ 
-(2uUZ^Mg  -  MJg)  +  2u-Zg(zJV-ZgB))^  -  2u'((W  -  B)Mg , 
+  (  x^W  -  XgB  )Zg  +  ((iV  -  B)Mg  +  (  XfjfV  -  XgBjZg  *  =0 


Solution  Set  (CD-0.0,  Zgb-0.1,  Xgb/L- -0.0001) 
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Figure  4-1.  Exact  Solution  Set  for  Xgb/L--0.0001. 


where  Q  now  equals  the  terms  -pQ  present  in  equation  4.3.  Defining  the  following 


parameters 

=  (4  8) 

Xc8=Xo-Xb<  (4.9) 

-G8=‘G-‘S.  (4  10) 

z^W-2,B  =  :^,fV-2,S,  =  &„  (4  11) 

XaiV-XgB  =  XagfV~Xg6gS8^,  (4.12) 

A.  =  l  Fn\  (4  13) 

and  substituting  them  into  equation  4  7  will  allow  the  coefficients  in  front  of  the  various 
powers  of  w  to  be  written  as  follows, 

U’  A,  =  Z,{2^glV-2g5gJ  =  Zg(\u'lV-2g8g),  (4  14) 

w|w|;  A^  =  2u'CaA(Mg-x^Zg),  (4.15) 

H'  A.=-2u'(y'(Z^Mg-MJ^g}  +  Zg(‘ku'W-2ghB)),  (4.16) 

4  =  -2u-(~hgMg -Zg(x^gW-Xghg))  (4. 17) 

W-:  A^=(-8gMg-Zg(XcgW-Xghg)).  (4  18) 

The  resulting  form  of  equation  4  7  can  be  written 

/44H'’  +  /43H'|w|  +  ,4;M'  +  /4| +  /4|H'*  =0,  (4  19) 


which  is  a  generic  pitchfork  equation  If  it  is  assumed  that  the  submarine  is  neutrally 
buoyant,  and  equation  4  19  is  divided  by  Z^z^gW ,  then  the  above  coefficients  can  be 

redefined  as, 

A,  =  \,  (4  20) 


2uCoA(Mg-x^Zg) 


_g _ 


(421) 
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g^OB  2i,W 


(4.22) 


4= - 2g, 


(4.23) 


g^GB  " 


Letting  the  critical  parameter  \  be  defined  as 


(4.24) 


(425) 


will  allow  equation  4  19  to  be  put  into  the  following  pitchfork  solution  form; 

+  yXm'Im'I  -  lu'(\-\-  )w  +  2^t/‘X  +  PXh’*  =  0, 


where. 


u‘ 


r.g(Z,MB-MJB) 

^  ^7  Tt/  ' 


(426) 


(4.27) 


(4  28) 


(429) 


To  determine  if  the  above  manipulation  and  definitions  are  correct  and  valid,  take  the 
symmetric  case  that  was  developed  in  the  previous  chapter,  i  e  ,  =0  This  causes 


equation  4  26  to  reduce  to 


yv(  w‘  +  yXIh-I  -  2«V1  +  ^X>>  =  0, 


(4  30) 


resulting  in  a  solution  always  occurring  at  m'=0,  with  two  more  solutions  appearing  at  the 


bifurcation  value 


i-kCX„  =  o 


(4.31) 


Rearranging  and  substituting  the  appropriate  values  into  equation  4  31  will  yield 


which  is  identical  to  equations  3  9  and  3.20,  thus  verifying  the  manipulation 

To  find  the  relationship  between  Fn  and  xq^  the  function  h(w,X,^),  and  its  partial 
derivative  must  be  determined  These  two  functions  must  then  be  set  equal  to 

zero  and  to  each  other  in  order  to  obtain  a  function  independent  of  w.  As  seen,  equation 
4  30  is  much  too  complex  to  complete  these  requirements,  but  by  neglecting  the  terms 
Xh'|h'|  and  Xm  '  in  equation  4.26  a  simplified  pitchfork  of  the  form, 

w'-2MVl  +  ;X>  +  2(itt'X  =  0,  (433) 

or  in  general  terms  h(w.X.^J  =  0,  can  be  obtained  which  can  be  manipulated  to  meet  the 
above  mentioned  requirements,  thereby  determining  the  critical  (X,^J  curve  By  talcing 
the  panial  derivative,  of  equation  4.33,  the  following  equation  results, 

A.  =3w=-2uVl  +  CX;  (4  34) 

This  equation  can  be  set  equal  to  zero,  yielding 

w-=y«Vl  +  ;XA  (4  35) 

which  can  be  substituted  into  equation  4  33  to  obtain  the  value  of  w  where  the  function 
and  it's  partial  derivative  is  equal  to  zero; 


3PX 

w  = - . 

2n+;x; 


(4  36) 


This  value  of  h  can  now  be  substituted  back  into  equation  4  35  to  obtain  the  desired 
f(X.^).  or 

27p-X-=8ttVl  +  ;X/  (4  37) 

Attempting  to  solve  this  equation  analytically  for  X  as  a  function  of  P ,  to  determine 
the  shape  of  this  curve  is  difficult  Therefore,  an  approximation  to  this  equation  is  needed 
for  small  values  of  p  By  rearranging  equation  4  37  into  the  form 


an  equation  that  determines  |3  as  a  function  of  X  is  now  produced  If  this  function  is 
represented  in  a  Taylor  series  expansion  and  simplified,  then  equation  4  38  can  be  written 
as 

=  (4.39) 

which  can  be  considered  as  an  approximation  to  the  (X.^J  curve,  with  equation  4  37  as 
the  exact  (X.^)  curve  Through  inspection  it  can  be  seen  that  the  general  shape  of  this 
curve  is  that  of  a  cusp 

Figure  4-2  shows  a  comparison  of  four  bifurcation  curves  The  curves  represent  the 
relationship  between  the  critical  speed  f/,  which  can  be  converted  into  the  non  dimensional 
parameter  X  using  equation  4  25,  and  the  value  of  xq^,  which  also  can  be  converted  into  a 
non  dimensional  parameter  p  using  equation  4.27  The  curves  labeled  exact  bifurcation 
set,  pitchfork  bifurcation  set,  exact  cusp  and  approximate  cusp  come  from  the  numerical 
solution  to  equations  4  7,  4.26,  4.37  and  4  39  respectively  Observe  that  in  the  general 
vicinity  of  XQg=0  and  near  the  critical  speed  the  shape  of  each  of  these  curves  is  indeed 
that  of  a  cusp,  with  the  peak  occurring  at  the  critical  speed  of  1  9  fps  Now,  with  the 
comparison  conducted  in  Figure  4-2  complete,  equation  4  7,  the  exact  pitchfork  solution 
will  be  used  for  the  remainder  of  this  analysis 

Returning  to  equation  4  7,  it  can  be  seen  that  there  is  only  one  additional  parameter, 
Cf),  that  the  cusp  curve  produced  using  equation  4  39  does  not  incorporate  (because 
and  XQg  are  incorporated  into  the  non  dimensional  parameters  X  and  P)  It  is  therefore 
necessary  to  conduct  a  sensitivity  analysis  to  determine  the  effect  that  Cq  has  on  the  shape 
of  the  cusp  curve  Figure  4-3  displays  the  results  obtained  when  this  sensitivity  analysis 
was  conducted 
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As  seen  in  the  plot,  the  drag  coefficient  has  no  effect  on  the  location  of  the  peak  of 
the  cusp  The  drag  coefficient,  however,  does  modify  the  slope  of  the  cusp  as  displayed  by 
the  four  cases  shown  in  Figure  4-3 

C  SOLUTION  SETS 

In  this  section,  the  equilibrium  solutions  to  equations  4  3  and  4  5  will  be  investigated 
As  determined  in  the  previous  section,  the  critical  parameters  that  control  the  location  and 
shape  of  the  pitchfork  solution  sets  are  X ,  P  and  the  drag  coefficient  Q  By  running  the 
FORTRAN  program  of  Appendix  A  for  several  values  of  L  the  shape  and  trend  of 
the  solution  sets  can  be  determined  The  result  of  this  comparison,  for  a  Cd  =  0.0,  is 
displayed  in  Figure  4-4  This  plot  supports  the  symmetry  about  x^g  =0  displayed  in  the 
earlier  cusp  curves  Observe  that  the  critical  or  bifurcation  point  moves  to  a  lower  value  of 
X  as  xqq  increases,  and  that  for  low  values  of  X  the  equilibrium  solutions  converge  to  the 
same  value  independent  of  the  value  of  xq^. 

When  the  comparison  is  made  regarding  the  effect  that  the  drag  coefficient  has  on  the 
solution  sets,  for  a  given  value  of  XQg,  Figure  4-5  is  produced  This  variation  in  drag 
coefficient  for  a  constant  value  of  Xq^  has  a  similar  effect  on  the  movement  of  the 
bifurcation  point  as  fixing  the  drag  coefficient  and  varying  xq^  The  difference  is,  however, 
evident  by  studying  the  stable  solutions  The  stable  solutions  for  non-zero  drag 
coefficients  become  more  linear  as  the  value  of  increases  This  results  from  the 
increased  effect  of  the  w'|u'|  term  in  equation  4  3  as  the  value  of  Cp  increases 

Although  Figure  4-5  shows  the  actual  solution  to  the  exact  pitchfork  equation  it  is 
not  realistic  since  the  effect  of  dive  plane  saturation  is  not  considered  Once  the  dive 
planes  saturate,  the  equations  used  to  obtain  the  steady  state  solutions  displayed  in 
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Figure  4-4.  Comparison  of  Exact  Solution  Set  for  Different  Values  of  Xgb/L. 
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Figures  4-4  and  4-5  are  no  longer  valid  Instead,  equations  4.1  and  4.2  are  modified,  still 
remembering  that  the  assumption  of  neutral  buoyancy  remains  in  affect,  resulting  in  the 
following  two  equations, 

-  Q^h'Ih’I  +  W-Zg5^,  =  0,  (4.40) 

Mjuw-C,^x^Av\yii\-WxQgCosQ-ZQgW  sin%-¥  =0.  (4.41) 

To  obtain  an  equation  that  can  be  used  to  compute  the  steady  state  6  solution  while 
taking  into  consideration  dive  plane  saturation  requires  solving  equation  4.40  for  w  and 
then  substituting  this  value  into  equation  4  41  which  should  yield  at'  equation  independent 
of  H  The  ability  to  perform  this  manipulation  is  hindered  by  the  inclusion  of  the  absolute 
value  term  in  both  of  the  above  equations  However,  if  the  absolute  value  terms  are 
neglected  by  considering  the  case  of  C[y=0.Q,  then  a  single  equation  for  the  steady  state 
value  of  0  may  be  obtained 

Neglecting  the  drag  affects,  and  performing  the  algebra  will  yield  an  equation  of  the 
following  form, 

( A/.Zj  -  A/gZJ«-5„,  +  Z„irrXogCO50  +  Zoa5me;  =  0,  (4  42) 

which  may  be  used  to  determine  the  values  of  6  once  the  dive  planes  have  saturated.  The 
validity  of  neglecting  the  drag  terms  is  questionable,  and  an  analysis  that  incorporates 
will  be  conducted  later  in  this  section  once  the  effect  of  varying  Xq^  is  determined  Figure 
4-6  displays  the  changes  that  occur  to  the  solution  set  of  Figure  4-1  when  saturation  is 
considered  Notice  that  saturation  occurs  when  the  solution  set  is  at  approximately  nine 
degrees,  and  that  when  multiple  solutions  occur,  the  stable  equilibrium  states  are 
represented  by  the  saturation  portions  of  the  solution  set 
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Figure  4-6.  Exact  Solution  Set  for  Xgb/Lr=-0.0001  With  Saturation  Included. 
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When  the  saturation  condition  is  applied  to  the  previously  displayed  solution  sets  for 
different  values  of  xq^  (Figure  4-4),  the  results  shown  in  Figure  4-7  are  produced.  Items 
of  interest  in  Figure  4-7  are  that: 

1 )  The  bifurcation  point  on  each  solution  set  moves  to  a  lower  value  of  X  as 
the  absolute  value  of  increases 

2)  The  range  where  multiple  solutions  exist  decreases  as  the  absolute  value  of 
XQg  increases 

3)  The  amount  of  trim  produced,  in  still  water,  by  the  given  value  of  XQg  can  be 
obtained  from  the  point  where  the  upper  an  lower  solution  branches 
converge 

Incorporating  the  saturation  condition  into  the  solution  sets  raises  some  additional 
concerns  To  begin,  wiiat  if  the  point  of  dive  plane  saturation  occurs  after  the  bifurcation 
points  identified  in  Figure  4-4“^  If  this  condition  did  occur,  referring  to  Figure  4- 1  which 
displays  the  stable  and  unstable  solution  branches,  there  would  be  at  least  one  value  of  X 
where  there  were  three  stable  solution  branches.  To  determine  if  this  condition  is  possible 
a  comparison  of  solution  sets  with  and  without  saturation  is  necessary.  This  comparison  is 
shown  in  Figure  4-8  As  this  plot  depicts,  the  range  between  the  bifurcation  point  and 
saturation  point  increases  as  the  value  of  xqb  increases  This  result  will  allow  the  fore 
mentioned  concern,  of  three  stable  solution  branches,  to  be  disregarded 

A  second  concern  is  the  validity  of  the  previously  developed  bifurcation  cusp  curves 
Since  the  exact  bifurcation  equation,  equation  4-7,  does  not  include  saturation,  the  cusp 
curves  produced  using  this  equation  will  not  accurately  predict  the  occurrence  of  multiple 
solutions  However,  if  equation  4  42  is  expanded  in  a  Taylor  series,  it  can  be  manipulated 
into  a  form  that  will  allow  a  cusp  curve,  which  includes  saturation,  to  be  developed 
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Figure  4-8.  Comparison  of  Bifurcation  Points  and  Dive  Plane  Saturation  Points. 


Performing  this  manipulation  will  produce  the  following  equation; 


(4,43) 


(Z,M^-M,zj(hj\HZBZJK.r) 

By  converting  u-  and  xq^  into  the  non  dimensional  parameters  X  and  ^  the  cusp 
curve  shown  in  Figure  4-9  can  be  produced.  This  cusp  displays  the  saturation  point  for 
both  solution  branches,  and  as  seen  it  differs  quite  extensively  from  the  exact  bifurcation 
cusp  without  saturation.  If  Figure  4-9  is  replotted  to  allow  the  peak  area  of  the  cusp  to  be 
more  accurately  represented,  then  Figure  4-10  is  produced.  By  using  these  two  cusp 
curves  then  the  general  form  of  the  resulting  solution  set  for  any  path  through  the  cusp  can 
be  predicted  It  is  this  prediction  that  the  subsequent  section  will  address 


D.  PATH  FORMULATION 

The  ability  to  accurately  predict  how  the  vessel  will  respond  to  various  changes  in 
operating  conditions  is  critical  for  proper  control  of  the  vessel  It  is  this  prediction  that 
necessitates  the  need  for  path  analysis  to  be  conducted  If  the  ballast  condition  on  the 
vessel  is  held  constant,  Xq^  =-0.0001,  and  the  speed  of  the  vessel  is  varied,  path  A  of 
Figure  4-11  is  obtained  Neglecting  for  the  moment  the  saturation  cusp  shown  in 
Figure  4-11,  it  would  be  expected  that  a  single  steady  state  solution  would  exist  for  all 
values,  of  the  parameter  sqrt(X),  greater  than  the  point  where  the  path  intersects  the  solid 
cusp  curve  At  the  speed  where  the  path  intersects  the  cusp  the  bifurcation  point  occurs, 
and  for  speed  values  which  place  the  path  inside  the  cusp  multiple  steady  state  solutions 
occur  The  results  of  this  path  analysis  are  depicted  in  Figure  4-1,  presented  earlier  in  this 
chapter 

If  attention  is  returned  to  the  path  through  the  saturation  cusp  of  Figure  4-11,  then 
the  same  type  of  predictions  made  for  the  path  through  the  solid  cusp  can  be  made  for  the 
this  cusp,  with  some  slight  differences  As  the  speed  decreases  the  path  intersects  the 
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upper  portion  of  the  dotted  cusp  at  an  approximate  value  of  1.08  At  this  value  saturation 
occurs  on  one  of  the  solution  branches  but  there  remains  only  one  steady  state  solution 
since  the  path  has  not  yet  entered  the  cusp  As  the  speed  continues  to  decrease,  the  path 
intersects  and  enters  the  saturation  cusp  at  an  approximate  value  of  0.99  with  saturation  of 
the  other  solution  branch  and  the  existence  of  multiple  steady  state  solutions  occurring. 
These  results  are  evident  in  Figure  4-6 

Now,  if  the  speed  of  the  vessel  is  held  constant,  and  the  loading  conditions  on  the  hull 
are  changed,  then  the  paths  shown  in  Figure  4-12  are  obtained  Path  B,  at  a  Froude 
Number  of  1  1 ,  is  outside  the  cusp,  therefore  only  a  single  steady  state  solution  for  each 
value  of  XQg  L  should  exist  While  path  C,  at  a  Froude  Number  of  1  0,  passes  through  the 
cusp  at  ±  0  0001  values  of  xq^L  resulting  in  single  steady  state  solutions  for  values  less 
than  or  greater  than  ±  0  0001,  while  values  between  ±  0  0001  result  in  multiple  steady 
state  solutions  Through  numerical  computation  the  solution  set  as  a  function  of  XQg  L  is 
obtained  with  the  results  shown  in  Figure  4-13 

As  seen  in  Figure  4-13,  path  B,  represented  by  the  dotted  line,  does  in  fact  result  in  a 
single  steady  state  solution  throughout  the  range  of  Xq^  /,.  however,  for  path  C  this  is  not 
the  case  Path  C.  represented  by  the  solid  line,  exhibits  multiple  solutions  between  the 
XQgL  values  of  ±  0  0001,  as  predicted  This  is  a  hysteresis  curve,  with  the  unstable 
solutions  occurring  between  the  comers  of  the  curve  If  the  value  of  xq^L  began 
at  -0  0005  and  progressed  to  the  positive  value  of  0  0001,  the  solution  set  would  remain 
on  the  curve  with  0  taking  on  the  values  shown  Once  the  value  of  XQg  L  becomes  greater 
than  0  0001  the  solution  jumps  from  -9°  to  7  5°  This  same  jump  would  be  evident  if  the 
value  of  XQg  L  began  positive  and  progressed  negative  These  jumps  indicate  that  if  the 
dive  planes  are  held  constant,  then  the  vessel  will  instantaneously  change  its  pitch  angle 
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73 


To  counter  this  effect  the  dive  plane  angle  must  be  reversed  It  is  this  dive  plane  reversal 
that  the  following  chapter  will  discuss. 

Before  proceeding  to  the  next  chapter  it  is  necessary  to  determine  the  effect  that  drag 
has  on  the  previously  developed  bifurcation  cusps  Using  equations  4  40  and  4.41,  and 
substituting  equation  4  4  into  both  equations  will  result  in  the  following  two  equations, 

Z,^«^tan0-Qy4«’ tan6|tan6|  +  ZjM^5„,  =  0,  (4  44) 

M^uw-CpX^Au^  tan€lft&n€^-fVX(;gCOs6- Zf;gWsind+  =0  (4  45) 

By  using  the  MATLAB  function  fzeros  contained  in  the  program  of  Appendix  A,  the 
value  of  6  that  solves  equation  4  44  can  be  determined  as  a  function  of  speed  This 
relationship  can  then  be  used  to  solve  equation  4  45  to  obtain  the  relationship  between  Xfjg 
and  u 

With  this  relationship  obtained  the  cusp  curves  of  Figure  4-14  can  be  plotted  It  can 
be  seen  that  for  a  given  value  of  xq^L  increasing  the  drag  coefficient  forces  the 
bifurcation  point  to  a  lower  speed  An  interesting  item  to  note  is  the  diamond  shape  peak 
on  these  cusps  This  is  a  result  of  saturation  not  occurring,  for  a  small  range  of  speeds,  in 
the  small  region  around  X(~gL=0  Path  formulation  through  these  cusps  can  be  developed 
as  shown  earlier  in  this  section 
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V.  DIVE  PLANE  REVERSAL 


A.  STERN  PLANE  REVERSAL 

The  ability  of  a  submarine  to  avoid  detection  and  perform  its  mission  is  dependent  on 
its  ability  to  maintain  ordered  depth.  Modem  asymmetric  submarines  exhibit  an  inherent 
phenomenon  known  as  dive  plane  reversal,  which  is  difficult  to  deal  with.  To  introduce 
this  physical  occurrence,  refer  to  Figure  5-1  [Ref  14]  and  consider  the  case  of  wanting  to 
dive  the  vehicle  to  a  deeper  depth.  At  moderate  to  high  speeds  the  stem  planes  would  be 
deflected  by  an  angle  5^,  this  would  in  turn  cause  a  force  and  moment  to  be  developed  due 
to  the  angle  of  attack  of  the  planes.  The  vessel  will  respond  by  pitching  downward,  and 
will  assume  some  angle  of  attack  to  the  fluid  flow  causing  hull  forces  and  moments  to 
develop  The  stem  plane  and  hull  forces  bring  the  ship  to  some  steady  state  pitch  angle 
which  is  obtained  when  the  restoring  moment  and  pitching  moments  balance.  The  vessel 
now  has  the  ability  to  fly  itself  to  the  ordered  depth  with  the  normal  hull  force  assisting  in 
pushing  the  vessel  to  a  deeper  depth 

At  low  speed  we  would  expect  the  same  to  be  true,  however  it  is  not  If  the  stem 
planes  are  set  to  an  angle  5^ ,  this  would  cause  a  force  and  moment  to  be  developed  due  to 
the  angle  of  attack  of  the  planes,  the  vessel  will  pitch  downward  creating  a  pressure  field 
around  the  vessel  due  to  the  angle  of  attack  with  the  fluid  resulting  in  hull  forces  and 
moments  to  develop  The  vessel  once  again  steadies  out  at  some  steady  state  pitch  angle 
when  the  restoring  moment  and  pitching  moments  balance  Since  the  speed  of  the  vessel  is 
not  as  great  as  in  the  previous  case,  the  hull  forces  and  moments  are  not  as  large,  therefore 


instead  of  the  vessel  fving  itself  to  the  ordered  depth  with  the  hull  force  assisting,  the 
stem  plane  force  actual  pushes  the  vessel  upward  causing  it  to  rise  instead  of  dive,  with 
some  non-zero  pitch  angle.  [Ref  14]  This  condition  was  evident  in  the  simulation  shown 
in  Figure  3.3 

With  this  description  of  the  vehicle  response  complete,  it  is  obvious  that  to  dive  the 
vessel  at  low  speeds,  the  opposite  stem  plane  angle  must  be  ordered.  This  will  result  in  the 
vessel  being  pushed  to  the  ordered  depth  at  some  bow  up  attitude  This  required  shift  in 
dive  plane  angle  is  what  is  known  as  STERN  PLANE  REVERSAL 

To  tie  together  the  previously  completed  bifurcation  analysis  with  the  physical 
occurrence  just  discussed.  Figure  5-2  is  produced  Figure  5-2  presents  the  pitch  angle  and 
stem  plane  angle  required  for  straight  and  level  flight  with  the  bow  planes  centered  As  the 
vessel's  speed  decreases,  the  metacentric  moment  ( MJd)  makes  it  increasingly  difficult  to 

maintain  a  positive  pitch  angle  on  the  hull  The  stem  plane  angle  is  negative  trying  to 
produce  a  pitch-up  moment,  but  the  downward  force  of  the  stem  planes  requires  an  even 
larger  pitch  angle  [Ref  14]  Between  the  speeds  of  I  08  and  I  18  knots,  the  required  stem 
plane  angle  is  beyond  the  normal  operating  range  and  the  vessel  cannot  be  flown  straight 
and  level  Below  the  cntical  speed,  the  stem  plane  and  vessel  pitch  angle  must  reverse  sign 
to  maintain  neutral  tnm 

In  a  final  effort  to  show  the  physical  significance  of  stem  plane  reversal  the  force 
input  to  the  stern  planes  is  plotted  If  the  control  gains  ate  allowed  to  change  as  the  speed 
of  the  vessel  changes  figure  '-t  is  produced  This  plot  sh<»ws  the  norr  dimensional  force 
per  unit  depth  error  input  tc'  the  stern  plane  as  a  tunciion  of  speed  li  is  basicalK  the 
control  system  gam  used  to  set  the  stern  plane  angle  4s  can  he  seen  as  the  critical  speed 
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is  approached  the  value  of  input  tends  to  infinity  requiring  the  stem  planes  to  operate 
outside  normal  values.  As  the  speed  passes  the  critical  value  the  sign  of  this  input  changes 
requiring  the  stem  planes  to  operate  in  the  opposite  direction  This  reversal  in  sign  is 
evident  of  the  stem  plane  reversal  phenomenon  shown  and  discussed  earlier  in  this 
chapter 

B.  BOW  PLANE  REVERSAL 

The  reversal  effect  that  was  shown  to  occur  with  the  stem  planes  can  also  be  shown 
to  occur  with  the  bow  planes.  If  the  value  of  a  is  allowed  to  tend  to  oo,  indicating  no  stem 
plane  control  and  only  bow  plane  control,  then  the  same  analysis  conducted  throughout 
this  thesis  for  stem  plane  reversal  can  be  performed  for  bow  plane  reversal.  If  the  critical 
Froude  number,  (sqrt(X)),  is  plotted  as  a  function  of  a,  then  the  results  shown  in 
Figure  5-4  are  produced  This  plot  shows  that  for  bow  plane  control  only,  the  critical 
value  approaches  2  05  as  a  approaches  oo 

To  show  that  the  same  type  of  results  can  be  obtained  for  the  bow  plane  analysis,  the 
case  of  neutral  buoyancy,  =  0  and  dive  plane  saturation,  with  bow  plane  control  only 
was  investigated  Referring  to  Figure  5-5,  it  can  be  seen  that  at  a  Froude  number  of  2  05 
multiple  steady  state  solutions  occur  These  results  are  similar  to  the  results  displayed  in 
Figure  -L5  The  major  differences  are  that  the  critical  point  occurs  at  a  Froude  number  of 
2  0^  vice  1  os  and  the  maximum  value  of  0  is  only  ^  4  5°  vice  ±9  5°  The  reason  for 
these  changes  is  due  to  the  values  of  and  being  only  one-half  the  respective  value 

for  the  stern  planes 
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STEADY  STATE  VALUES  OF  THETA  vs.  FROUDE  NUMBER 


C.  BIAS  EFFECTS 


As  a  final  area  of  analysis,  the  condition  of  operating  near  the  surface,  such  as  during 
periscope  operations,  must  be  investigated  to  determine  the  effects  this  will  have  on  the 
previously  developed  bifurcation  cusps  Operating  near  the  surface  will  require  the 
incorporation  of  a  suction  force  and  moment  into  the  equations  of  motion  This  force  and 
moment  is  a  function  of  the  distance  away  from  the  surface,  and  the  angle  the  vessel 
makes  with  the  surface.  Therefore,  as  a  simple  model  of  suction  effects  consider  the  case 
of  operating  with  excess  buoyancy  To  conduct  the  analysis  requires  the  use  of  equations 
4  1,  4  2  and  4  4  Performing  the  same  substitution  and  solution  technique  as  was  used  in 
Chapter  V  to  obtain  the  saturation  curves  that  incorporated  the  drag  coefficient  will  result 
in  solving  the  following  equations  to  obtain  the  bifurcation  cusps  that  would  quasi-model 
surface  effects, 

Z,uw-^f>CfyAu‘tSLn6\un6^-¥{iV-B)cosd-*-Zfu'5^^,  =  0.  (5  1) 

M.um  tan6(tan6j-(x,W'-x^^)cos0 

-  r^5)sin0+  =  0 

If  the  value  of  (  '/)  is  set  equal  to  zero  then  the  saturation  cusp  produced  in 
Figure  4-10  can  be  ompared  to  the  results  obtained  for  several  cases  of  excess  buos  ancs 
Since  It  IS  the  saturation  cusp  that  controls  the  location  of  the  bifurcation  point,  as  shoun 
in  previous  the  chapter  the  non  saturation  cusp  vsill  be  neglected  tor  the  remainder  of  the 
analvsis  Figure  shosss  the  results  obtained  as  excess  buovancs  is  increased  As  can  be 
seen  the  cusp  peak  mosrs  to  dovsn  and  U'  the  right  occurring  ai  a  loccer  speed  and  at  a 
non  zero  value  of  r  /  The  bods  I’t  the  cusp  als«' shifts  to  the  right  as  seen  vsith  the 
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entire  cusp  occurring  in  the  positive  region  of  the  range  for  a  value  of  o  o*.  ®o 

excess  buoyancy 

To  determine  how  the  value  of  drag  coefficient  effects  the  cusp  curves  consider  the 
case  of  Cp,  =  0  3  shown  in  Figure  5-7  These  cusps  were  produced  for  the  same  valuev  I't 
excess  buoyancy  used  in  Figure  5-6  Notice  that  the  cusp  curves  shift  lo  the  nghi  hv  ihe 
same  value  seen  in  Figure  5-6,  but  the  peak  value  now  occurs  at  an  even  a  lower  vpeed 
In  order  to  observe  the  sensitivity  of  the  biased  cusp  to  drag  onK  consider  the  i-a^t 
of  setting  the  excess  buoyancy  to  0  01  ®o  and  varvmg  the  drag  coefficient  shmsi 
Figure  5-8  These  results  indicate  that  increased  drag  causes  the  cntiLai  speed  tt-  i>ni'  a’ 
a  lower  value  They  also  show  that  drag  tends  to  counter  the  biasing  cfte«.t  tha'  fv>t' 
buoyancy  was  previously  show  to  have  This  is  evident  bv  the  peak  point  dnhinc  •  o,, 
left  as  the  value  of  drag  coefficient  increa.ses 

In  summary,  it  can  be  stated  that  operating  near  the  surtacc  will  lenc  •  edu. . 
critical  speed  determined  in  previous  chapters  and  the  vessel  will  have  i  rn  pia.  r 

different  tnm  condition  to  counter  surface  effects  and  that  incieaMnc  l  ac  a  ^  . . 

the  critical  speed  but  tends  to  counter  the  tnm  condition  leguned 
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Figure  5-7.  Bias  Effects  Caused  by  Operating  Near  the  Surface  Considering  Drag. 


sqrt(  Lambda) 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

•  A  method  for  analyzing  submarine  depth  keeping  and  dive  plane  reversal  has  been 
developed  and  presented.  This  analysis  identified  the  three  phenomena,  a  positive  real 
eigenvalue,  a  steady  state  pitchfork  bifurcation,  and  an  uncontrollable  system,  that  lead  to 
vessel  instability. 

•  The  critical  parameters  have  been  identified  that  control  the  bifurcation.  These 
parameters  combine  speed,  loading  conditions  and  hydrodynamic  coefficients  into 
non-dimensional  values  that  may  be  applied  to  other  models  or  vessels 

•  Since  the  non-dimensional  Froude  numbers  are  based  on  the  metacentric  height 
and  LCG-LCB  separation,  Froude  scaling  can  be  used  to  apply  this  analysis  to  full  scale 

vessels,  allowing  for  the  reduction  of  model  testing  during  design  phases. 

*  * 

*  ^ 

B.  RECOMMENDATIONS 

Some  recommendations  for  further  research  are  as  follows: 


•  Develop  and  incorporate  into  the  bifurcation  analysis  a  more  complete  expression 
for  the  free-surface  effects 

•  Conduct  a  Hopf  bifurcation  analysis  to  identify  the  critical  parameters  that  lead  to 
vehicle  instability  at  high  speeds 


90 


•  ••••••• 


APPENDIX  A 


COMPUTER  PROGRAMS 

%  THIS  PROGRAM  IS  THESIS  I  M  IT  SIMULATES  THE  RESPONSE  OF  THE 
DTRC  SUBOFF 

%  MODEL  IT  PLACES  THE  POLES,  CALCULATES  THE  CRITICAL  SPEED.  AND 
USES 

%  PROGRAM  THESISDE  M  TO  SIMULATE  THE  VERTICAL  PLANE  REPONSE 


global  zg  m  zgb  u  U  Zq  Zw  Zdlt  Mq  Mw  B 1  Mdlt  K1  M  bm  xL  Zval  Mval  rho 

%  ESTABLISH  GEOMETRIC  PARAMETERS 
W=1556  2363,B1=1556  2363; 

Iy=561.32, 
g=32.2; 
m=W/g, 
rho=l  94, 

L=  13  9792, 
xg=0,zg=0. 1 , 
zgb=0  1, 

%  BEAM  AND  LENGTH  DEFINITION  FOR  TRAPAZIODAL  RULE  USE 
bm=[0  485  658  778  871  .945  1  01  1  06  118  1  41  1.57  1.66  1.67  1  67  1  67  1  63 
1  37  919  448  195  188  168  132  053  0], 

xl=[0  1  2  3  4  5  6  7  1  2  3  4  77143  10  15  1429  16  17  18  19  20  20  1 
20  2  20.3  20  4  20  4167], 
xI=(L/20)*xl, 
xL=xl-L/2, 


%  FACTORS  TO  CONVERT  FROM  NON-DIMENSIONAL  VALUES  TO 

DIMENSIONAL  VALUES 

ndl=  5*rho*L^2, 

nd2=  5*rho*L^3, 

nd3=  5*rho*L^4; 

nd4=  5*rho*L'5, 


w 


Alpha=0; 

%  NON-DIMENSIONAL  HYDRODYNAMIC  COEFFICIENTS 
Zqdnd=-6 . 3  3  e*4;Z  wdnd=- 1.45 29e-2  ,Zqnd=7 . 545e-3  ;Zwnd=- 1  3  9 1  e-2 , 
Zds=(-5  603e-3);Zdb=0  5*(-5  603e-3);Zdltnd=(Zds+Alpha*Zdb); 
Mqdnd=-8  8e-4;Mwdnd=-5 . 6 1  e-4,Mqnd=-3  702e-3  ,Mwnd=  1  0324e-2. 
Mds=(-0.002409);Mdb=0.5*(0.002409),Mdltnd=(Mds+A!pha*Mdb), 

%  CONVERTION  OF  COEFFICENTS  INTO  DIMENSIONAL  VALUES 
Zqd=nd3  *Zqdnd,Zwd=nd2'''Zwdnd,Zq=nd2*Zqnd,Zw=nd  1  *Zwnd, 
Zdlt=ndl*Zdltnd, 

Mqd=nd4*Mqdnd;Mwd=nd3*Mwdnd;Mq=nd3*Mqnd,Mw=nd2*Mwnd, 

Mdlt=nd2*Mdltnd; 

%  FORMULATION  OF  A  AND  B  MATRIX  ELEMENTS 

Dv=(m-Zwd)*(Iy-MqdHm*xg+Zqd)*(ni*xg+Mwd), 

al  lDv=(Iy-Mqd)*Zw+(m*xg+Zqd)*Mw, 

al2Dv=(Iy-Mqd)*(m+Zq)+(m*xg+Zqd)*(Mq-m*xg), 

al3Dv=-(m*xg+Zqd)*W, 

bIDv=(Iy-Mqd)*Zdlt+(m*xg+Zqd)*Mdlt, 

a2 1  Dv=(m-Zwd)*Mw+(m*xg+Mwd)*Zw, 

a22Dv=(m-Zwd)*(Mq-m*xg)+(m*xg+Mwd)*(m+Zq). 

a23Dv=-(m-Zwd)*W, 

b2Dv=(m-Zwd)*Mdlt+(m*xg+Mwd)*Zdlt, 

al  I=al  lDv/Dv;al2=al2Dv/Dv.al3=aI3Dv/Dv, 
a2 1  =a2 1  Dv/Dv,a22=a22Dv/Dv,a23=a23Dv/Dv, 
b  1  =b  1  Dv/Dv,b2=b2Dv/Dv. 

%  ESTABLISHMENT  OF  POLE  LOCATIONS 
pl=[-  3  -  31  -  32  -  33], 


%  CALCULATION  OF  CRITICAL  SPEED 
Ucr=sqrt(-(a  1 3  *b2-a23  *b  1  )*zgb/(a  1 1  *b2-a2 1  *b  I )) 

%  INPUT  OF  SPEED  WHICH  THE  VESSEL  WILL  USE  IN  THE  SIMULATION 
U=input('enter  speed  of  vehicle  ') 

%  CALCULATION  OF  FROUDE  NUMBER 
Fn=sqrt(U^2/(zgb*g)) 


* 
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%  CALCULATION  OF  MATRICIES  FOR  POLE  PLACEMENT  AND 

CONTROLABILLITY  ANALYSIS 

u=5, 

a=[0  0  1  0,al3*zgb  al  l*u  al2*u  0;a23*zgb  a21*u  a22*u0,... 

-u  1  0  0], 

al=[0  0  1  0;al3*zgb  all*Ucr  al2*l4cr  0,a23*zgb  a21*Ucr  a22*Ucr  0. 
-Ucr  1  0  0], 


b=[0;bl*u^2;b2*u^2;0], 
b  1  =(0;b  1  •Ucr^2;b2*Ucr^2,0], 

Kl=place(a,b,pl); 

C=rank(ctrb(al,bl)), 

H=(Mw*Zdlt-Zw*Mdlt)/(zgb*B  I  *Zdlt); 
theta  1  =acos(H*U^2)*  1 80/pi, 

Dss=-(Zw/Zdlt)*tan(acos(H*U''2))*  1 80/pi; 

Z=(Dss+K  I  ( 1  )*acos(H*U^2)+K  U2)*U*tan(acos(H*U^2)))/ 

(-Kl(4)), 

%  ESTABISHES  VEHICLE  MASS  MATRIX 

M=[l  0  0  0,0  (m-Zwd)  -Zqd  0  ,0  -Mwd  (ly-Mqd)  0  ,0  0  0  1], 

%  SET  TIME  LIMITS  AND  INITIAL  CONDITIONS  FOR  ODE  SOLUTION 

t0=0,tfl=2500*L/U, 

x0=[0  0  0  1], 

%  SOLVES  VEHICLE  SIMULATION  ODE’S 
(t  1  ,X  I  ]=ode45('thesisde',t0,tfl  ,x0), 

%  REDEFINE  THE  OUTPUT  VECTOR  ELEMENTS  AND  FORM  DIVE  PLANE 
RESPONSE 

wl=Xl(  .2),ql=Xl(  ,3).zl=Xl(  ,4),thl=Xl(  .1), 
d=-(K  1  ( 1  )*th  1 +K 1  (2)*wl  +K 1  (3)*q  1 +K  l(4)*zl ). 

Il=length(d), 

for  n=)  11, 
if  d(n)>  4, 
d(n)=  4, 
elseif  d(n)<-  4. 

d(n)=-  4, 
end, 
end. 
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%  PLOT  THE  RESULTS 
subplot(21 1), 

plot(tl*U/L,thl*180/pi,tl*U/L,d*180/pi,'-');title(THET A/DELTA  vs  TIME  '); 
xlabelCnon-dimensional  time  ');y]abel('theta/delta  (degrees)');grid, 
gtext('— delta');gtext(['for  a  speed  of  num2str(U) '  fps']), 

plot(tl*U/L,zl/L);title('DEPTH  vs.  TIME');xlabei('non  dimensional  time'); 

ylabel('non  dimensional  depth  ');grid, 

pause, 

subplot(l  1 1), 

%  THIS  IS  PROGRAM  THESISDE  M  IT  CONTAINS  THE  SUBOFF  MODEL 
SYSTEM  OF 

%  EQUATIONS  TO  BE  SOLVED  BY  PROGRAM  THESIS  1  M 

%  SET  SYSTEM  OF  EQUATIONS  FUNCTION 
function  Xldot=thesisde(tl,Xl), 

%  SET  CONTROL  LAW 

dl=-Kl(l)*Xl(l)-Kl(2)*Xl(2)-Kl(3)*Xl(3)-Kl(4)*Xl(4), 

%  APPLY  DIVE  PLANE  SATURATION 
if  dl>  4, 
dl=4, 

eiseif  dl<-  4. 

dl=-  4, 
end. 

%  INCLUDE  DRAG  TERMS  AND  CALCULATE  AREA  AND  CENTROID 
CD=0, 

ifXl(2)==0, 

Xl(2)=le-5, 

end; 
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ifXl(3)==0, 

XU3)=le-5, 

end, 

Zval=bm  *((Xl(2)-xL  *Xl(3))/'3)  /(abs(Xl(2)-xL  '*X1(3))), 
MvaNbm  *(((X  I  (2)-xL  *X  I  (3 ))/3 )  *xL)  /(abs(X  I  (2)-xL  *X  I  (3 ))), 


ZDragval=0,MDragval=0, 
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%  trapaziodal  integration 
for  n=l  :length(xL)-l, 


Zclragval=0. 5  *(Zval(n)+Zval(n+  l))*(xL(n+ 1  )-xL(n)); 

Mdragval=0. 5  *(Mval(n)+MvaI(n+ 1  ))*(xL(n+ 1  )-xL(n)); 

ZDragvaNZDragval+Zdragval ; 

MDragval=MDragval+Mdragval; 

end, 

Zdr=(0  5)*rho*CD*ZDragval; 

Mdr=(0  5)*rho*CD*MDragval; 

%  VEHICLE  SYSTEM  OF  EQUATIONS 

F1=[X1(3),  m*(U*Xl(3)+zg*Xl(3r2)+Zq*U*Xl(3)+Zw*U*Xl(2)+U^2*ZdIt*dl-Zdr, 
m*(-zg*Xl(2)*Xl(3))+Mq*U*Xl(3)+Mw*U*Xl(2)- 
B 1  *zgb*sin(X  1  ( 1  ))+U^2*Mdit*d  1  -Mdr, 

-U*sin(X  1  ( I  ))+X  1  (2)*cos(X  1(1))]; 

mpr=inv(M), 

Xldot=mpr*FI, 


C  PRCKjRAM  steady 

c 

C  STEADY  STATE  SOLUTIONS  OF  SUBMARINES 
C  IN  THE  DIVE  PLANE  AT  LOW  SPEEDS 
C 

C  INPUT  DATA 
C 

C  ICON  =  1  COMPUTATION  OF  SOLUTION  SETS  (S  S) 

C  2  COMPUTATION  OF  BIFURCATION  GRAPHS  (B  G  ) 

C 

C  IV AR  =  1  Fn  VARIATION 
C  2  Xgb  VARIATION 

C 

REAL  MW,LENGTH,MASS,MDS.MDB.MDXAMBDA 
C 

DIMENSION  B(25).X(25),BX(25) 

C 

C  GEOMETRIC  PROPERTIES  AND  HYDRODYNAMIC  COEFFICIENTS 
C 

RHO  =  I  94 
LENGTH^  1 3  9792 
WRITE  (••)' ENTER  CDZ' 

READ  (•  •)CDZ 
CDZ  =CDZ*0  5*RHO 
BLO=  1556  2363 

WRITE! •.*)  ENTER  VALUE  OF  EXCESS  BUO' 

RE  AD  (*.*)  FAC 
WEIGHT=FAC*BUO 
MASS  =  WEIGHT/32  2 
XB  =00 
ZB  =00 

WRITE  (*,*)  ENTER  THE  MET  ACENTRIC  HEIGHT  ZGB' 
READ(*.*)ZG 

ZW  =-0  013910*0  5*RHO*LENGTH**2 
ZDS  =-0  005603*0  5*RHO*LENGTH**2 
ZDB  =-0  005603*0  25*RHO*LENGTH**2 
MW  =0  010324*0  5*RHO*LENGTH**3 
MDS  =-0  002409*0  5*RHO*LENGTH* *3 
MDB  =  0  002409*0  25*RHO*LENGTH**3 
C 

WRITE  (*,*)  'ENTER  THE  VALUE  OF  ALPHA 
READ  (*,*)  ALPHA 
ZD=ZDS+ALPHA*ZDB 


♦ 


•  • 


•  • 
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MD=MDS+ALPHA*MDB 


C 

C  DEFINE  THE  LENGTH  X  AND  BREADTH  B  TERMS  FOR  THE 

INTEGRATION 

C 

X(  1)=20.4167 
X(  2)=20.4000 
X(  3)=20.3 
X(  4)=20.2 
X(  5)=20  1 
X(  6)=20  0 
X(  7)=I9  0 
X(  8)=18  0 
X(  9)=17.0 
X(10)=16  0 
X(ll)=15  1429 
X(12)=10  0 
X(13)=7  7143 
X(14)=4  0 
X(15)=3  0 
X(16)=2  0 
X(17)=l  0 
X(18)=0  7 
X(19)=0  6 
X(20)=0.5 
X(21)=0  4 
X(22)=0  3 
X(23)=0.2 
X(24)=0  1 
X(25)=0  0 

C 

B(  1  >=0  00000 
B(  2)=0  03178 
B(  3)=0  07920 
B(  4)=0  10074 
B(  5)=0  1 1243 
B(  6)=0  11724 
B(  7)=0  26835 
B(  8)=0  55025 
B(  9)=0  81910 
B{10)=0  97598 
B(I  l)=l  00000 
B(12)=l  00000 
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B(13)=l  00000 
B(14)=0  99282 
B(15)=0  94066 
B(16)=0  84713 
B(17)=0  70744 
B(18)=0  63514 
B(19)=0  60352 
B(20)=0. 56627 
B(21)=0  52147 
B(22)=0  46600 
B(23)=0  39396 
B(24)=0  29058 
B(25  >=0  00000 
C 

DO  1  1=1,25 

B(1)=B(I)*  1  6667 
X(I)=(-X(I)+10  0)*LENGTH/20  0 
BX(I)=B(I)*X(1) 

I  CONTrWE 

CALL  TRAP(  1 8,B,X, AREA) 

CALL  TRAP(  1 8,BX,X,CH) 

XA=CHyAREA 

C  WRITE!*,*)  AREA  EQUALS '.AREA 
C  WRITE!*,*)  XA  EQUALS ',XA 
IFC=0 
C 

C  OPEN  RESULTS  FILES 

OPEN  !  10,FILE='R10  RES',STATUS='NEW') 
OPEN  !  1 1  ,FILE='R1 1  RES',STATUS='NEW) 
OPEN  !12,FILE=’R12  RES',STATUS='NEW') 
OPEN  !13,FILE='R13  RES’,STATUS='NEW') 
OPEN  !14,FILE=’R14  RES',STATUS='NEW') 
OPEN  !15,FILE='R15  RES',STATUS='NEW') 
OPEN  !16,FILE='R16  RES'.STATUS='NEW') 
OPEN  !20,FILE='C20  RES',STATUS='NEW’) 
OPEN  !21,FILE='S2I  RES'.STATUS='NEW  ) 
OPEN  !22,FILE='S22  RES'.STATUS='NEW  ) 
OPEN  !23,F1LE='S23  RES',STATUS=  NEW  ) 
OPEN  !:4.FILE='S24  RES',STATUS='NEW') 
OPEN  !31.FILE=’S3I  RES',STATUS='NEW') 
OPEN  !32,FILE='S32  RES'.STATUS='NEW') 
OPEN  !33,FILE='S33  RES'.STATUS='NEW') 
OPEN  !34,FILE='S34  RES',STATUS=’NEW  ) 
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OPEN  (4 1  ,FILE='C4 1  RES', STATUS=’NEW’) 
OPEN  (42,FILE='C42  RES’,STATL’S='NEW’) 
OPEN  (43,Fn/h>'C43  RES'.STATUS='NEW') 


READ  DATA  INTERACTIVELY 

WRITE  (MOOl)  • 

READ  (*,■*)  ICON 
IF  (ICON  EQ  2)  GO  TO  331 
WRITE  (M009) 

READ  (*,■»)  IVAR 

C}OTO(331,332),  IVAR  * 

331  WRITE  (M 002) 

READ  (*,*)  UMIN 
WRITE  (M003) 

READ  (*,*)  LTV1AX 

WRITE  (M 004)  ^ 

READ  (*.*)  ITER 

JL=ITER 

IF  (ICON  EQ  2)  GO  TO  332 
WRITE  (MOlO) 

READ  (*.*)  XG 
XG=XG*LENGTH 
GO  TO  300 

332  WRITE  (M005) 

READ  (*,'‘)  XGMIN 
UTGTE(M006) 

READ  (*,■•)  XGMAX  * 

VVTGTE  (M004) 

READ  (*,■*)  ITER 
JXG=ITER 

XGMIN=XGMIN*LENGTH 

XGMAX=XGMAX*LENGTH  • 

IF  (ICON  EQ  2)  GO  TO  334 
WRITE  (M012) 

READ  (*,*)  U 
334  CONTINUE 

(30  TO  300  • 

300  WRITE(M013) 

READ  (*,*)  ISET 
IF  (ICON  EQ  2)  (30  TO  370 
GO  TO  (350,360,400),  ISET 

• 


•  • 


•  • 
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c 

C  EXACT  SOLUTION  SET 
C 

350  DO  21=  LITER 

WRITE  (*.2000)  LITER 
GOTO  (351,352),  IVAR 

351  U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 

GO  TO  355 

3 52  XG=XGMIN+(XGMAX-XGMIN)*(I- 1  )/(ITER- 1 ) 

355  A1=ZW*MD-MW*ZD 

A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG-XB)*WEIGHT*ZD 
A4=  (ZG-ZB)*WE1GHT*ZD 
CALL 

EXS0LS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4,CDZ,AREA,ZW,ZD,ZG,MW. 
&  MD,BUO,XA) 

IF  (LL  GT  1)IFC=1 

2  CONTINUE 
GO  TO  5000 

C 

C  PITCHFORK  SOLUTION  SET 
C 

360  DO  3  1=1, ITER 

WRITE  (*,2000)  liter 
GO  TO  (361,362),  IVAR 

361  U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 

GO  TO  365 

362  XG=XGMIN+(XGMAX-XGMIN)*(I-1  )/(ITER-l ) 

365  A1=ZD*(ZG-ZB)*WEIGHT 

A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 
A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WE1GHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
A5=-ZD*(XG-XB)*WEIGHT 

CALL  PITSLS(IVAR,LL,IFC,U,XG,A1,A2,A3,A4.A5) 

IF  (LL  GT  1)IFC=1 

3  CONTINUE 
GO  TO  5000 

C 

C  SIMPLE  PITCHFORK  SOLUTION  SET 
C 

400  DO  401  1=1, ITER 

WRITE  (*,2000)  LITER 
GO  TO  (402,403),  IVAR 


•  • 


•  •  •  •  • 
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402  U=UMIN+(UMAX-UMIN)*(I-1)/(ITER-1) 

GO  TO  405 

403  XG=XGMIN+(XGMAX-XGMIN)*(I-1V(ITER-I) 

405  A1=ZD*(ZG-ZB)*WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
CALL  SIMSLS(IVAR,LL,IFC,U,XG,Al,A3.A4) 

IF(LLGT  1)IFC=I 
401  CONTINUE 
GOTO  5000 

370  GO  TO  (380.390.410.430),  ISET 
C 

C  EXACT  BIFURCATION  SET 
C 

380  DO  4  1=1. JXG 

WRITE  (*,2000)  l.JXG 

XG=XGMIN^(XGMAX-XGMIN)*(1- 1  )/(JXG- 1 ) 

IHS=0 
D0  5  J=1,JU 

U=UMIN+(UMAX-UMIN)*(  J- 1 )/( JU- 1 ) 

Al=  ZW*MD-MW*ZD 
A2=-CDZ*.AREA*(MD-XA*ZD) 
A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4=  (ZG*WEIGHT-ZB*BUO)*ZD 
CALL  EXN'UMS(NLrM.Al,A2,A3.A4,U) 

IF(JNE  1)  <30  TO  381 
UO=U 

NLfMO=NUM 
GOTO  5 

381  UN=U 
NUMN=NUM 

NDlF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  382 
UO=LT^ 

NUTV10=NUMN 
GOTO  5 

382  Ul=UO 
U2=UN 

DO  61  JJ=1.JU 

U=U1+(U2-U1)*(JJ-1)/(JU-1) 

Al=  ZW*MD-MW*ZD 

A2=-CDZ*AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 


Li 


•  • 


•  •  •  • 
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A4=  (ZG*WEIGHT-ZB*BUO)*ZD 
CALL  EXNUMS(NUM,A1.A2,A3,A4,U) 

IF  (JJ  NE  I)  GO  TO  383 
UO=U 

NUMO=NUM 
GO  TO  61 

383  LJN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  384 

UO=UN 

NUMO=NUMN 

61  CONTINUE 

GO  TO  5 

384  Ul=UO 
U2=UN 

DO  62  JJ=1,JU 
U=U  1  +(U2-U  1  1  )/(JU- 1 ) 

A1=ZW*MD-MW*ZD 
A2=-CDZ*  AREA*(MD-XA*ZD) 

A3=-(XG*WEIGHT-XB*BUO)*ZD+(WEIGHT-BUO)*MD 
A4=  (ZG*WEIGHT-ZB*BUO)*ZD 
CALL  EXNUMS(NUM,A1,A2,A3.A4,U) 

IF(JJNE  1)G0T0  386 
UO=U 

NUMO=NUM 
GOTO  62 
386  UN=U 

NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  385 

UO=UN 

NUMO=NUMN 

62  CONTINUE 

GOTO  5 

385  UPKUO+UN)/2  0 
IHS-IHS+1 

UPL=SQRT(UP/(32  2*LENGTH)) 

XGL=XG/LENGTH 

IPF=I0+IHS 

WRITE  (IFF,  100)  XGL.UP 

UO=UN 

NUMO=NUMN 


% 


4r 


•  • 
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L 


•  •  • 


•  • 


•  • 


•  • 


o  o  o 


5  CONTINUE 
4  CONTINUE 
GO  TO  5000 

PITCHFORK  BIFURCATION  SET 

390  DO  6  1=1, JXG 

WRITE  (*,2000)  I,JXG 

XG=XGMIN+(XGMAX-XGMIN)*(I- 1  )/(JXG- 1 ) 

IHS=0 
D0  7  J=1,JU 

U=UMIN+(UMAX-UMIN)*(J- 1 )/( JU- 1 ) 

Al=  ZD*(ZG-ZB)*WEIGHT 

A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 

A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 

AS=-ZD*(XG-XB)*WEIGHT 

CALL  PINUMS(NUM,A1,A2,A3,A4,A5,U) 

IF(J  NE  l)GOT0  391 
UO=U 

NUMO=NUM 

GOTO? 

391  UN=U 
NUMN=NUM 

NDIF=lABS(NUMO-NUMN) 

IF  CNDIF  EQ  2)  GO  TO  392 
UO=UN 
NUMO=NUMN 
GOTO? 

392  Ul=UO 
U2=UN 
DO?l  JJ=1,JU 

U=U  1  +(U2-U  1  )*(JJ- 1  )/(JU- 1 ) 

Al=  ZD*(ZG-ZB)*  WEIGHT 

A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 

A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 

A5=-ZD*(XG-XB)*WEIGHT 

CALL  PINUMS(NUM,AI,A2,A3,A4,A5,U) 

IF(JJ  NE  l)GOT0  393 
UO=U 

NUMO=NUM 

GOTO?l 


393  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  394 

UO=UN 

NUMO=NUMN 

71  CONTINUE 

STOP  1111 

394  Ul=UO 
U2=UN 

DO  72  JJ-1,JU 
U=U  1  +(U2-U  1  1  )/(JU- 1 ) 

Al  =  ZD*(ZG-ZB)*  WEIGHT 

A2=  2  0*U*U*U*CDZ*AREA*(MD-XA*ZD) 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
A5=-ZD*(XG-XB)*WEIGHT 
CALL  PINUMS(NUM,A1,A2,A3,A4,A5,U) 

IF  (JJ  NE  1)  GO  TO  396 
UO=U 

NUMO=NUM 
GOTO  72 
396  UN=U 

NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  395 

UO=UN 

NUMO=NUMN 

72  COmTNUE 

STOP  1112 

395  UP=(UO+UN)/2  0 
IHS=IHS+1 

UPL=SQRT(UP/(32  2*LENGTH)) 

XGL=XG/LENGTH 

IPF=10+IHS 

WRITE  (IPF.lOO)  XGL,UP 
UO=UN 
NUMO=NUMN 
7  CONTINUE 
6  CONTINUE 
GO  TO  5000 
C 

C  SIMPLE  PITCHFORK  BIFURCATION  SET 


410  DO  411  1=1,JXG 

WRITE  (*.2000)  1.JXG 

XG=XGMIN+(XGMAX-XGMIN)*(M  )/(JXG- 1 ) 

IHS=0 

00  412  J=1,JU 

U=UMIN+(UMAX-UMIN)*(  J- 1 )/( JU- 1 ) 

Al=  ZD*(2G-ZB)*WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
CALL  SINUMS(NUM,A1,A3,A4,U) 

IF  (J  NE  l)GOT0  413 
UO=U 

NUMO=NUM 
GO  TO  412 

413  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  414 
UO=UN 
NUMO=NUMN 
GO  TO  412 

414  Ul=UO 
U2=LrN 

D0  4I5  JJ=I.JU 
U=U1+(U2-U1)*(JJ-1)/(JU-1) 

Al=  ZD*(ZG-ZB)*  WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WE1GHT) 
A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
CALL  SINUMS(NUM,A1,A3,A4,U) 

IF  (JJ  NE  l)GOT0  416 
UO=U 

NUMO=NUM 

GOT0415 

416  UN=U 
NUMN=NUM 

NDIF=IABS(NUMO-NUMN) 

IF  (NDIF  EQ  2)  GO  TO  417 

UO=UN 

NUMO=NUMN 

415  -  CONTINUE 

STOP  1 1 1 1 

417  Ul=UO 


U2=UN 

D0418  JJ=1.JU 

U=U1+(U2-U1)*(JJ-1)/(JU-I)  • 

Al=  ZD*(ZG-ZB)*WEIGHT 

A3=-2*U*U*(U*U*(ZW*MD-MW*ZD)+ZD*(ZG-ZB)*WEIGHT) 

A4=  2  0*U*U*ZD*(XG-XB)*WEIGHT 
CALL  SINUMS(NUM,ALA3,A4.U) 

IF(JJ  NE  l)GOT0  419  • 

UO=U 

NLnvio=NLrM 

GO  TO  418 

419  UN=U 

NUMN=NUM  • 

NDIF=IABS(NUMO-NLTVfN) 

IF  (NDIF  EQ  2)  GO  TO  420 

UO=UN 

NUMO=NUMN 

418  CONTINUE  ^ 

STOP  1112 

420  UP=(UO+UN)/2  0 
IHS=IHS+1 

UPL=SQRT(UP/(32  2*LENGTH)) 

XGL=XG/LENGTH  . 

1PF=10+IHS 

WRITE  (IPF,  100)  XGL.UP 
UO=UN 
NUMO=NUMN 
412  CONTINUE 

411  CONTINUE  * 

GO  TO  5000 
C 

C  ANALYTIC  BIFURCATION  SET 
C 

430  WRITE  (M 01 7)  * 

READ  (*.•)  ICUSP 
GOTO  (432,433),  ICUSP 
C 

C  EXACT"  CUSP 

C  • 

432  DO  431  1=1, JL 

WRITE  (*,2000)  1,JU 
U=UMIN+(UMAX-UMIN)*(1- 1  )/(JU- 1 ) 

DELTA=32  2*(ZW*MD-MW*ZD)/(ZD*WEIGHT) 


•  •  • 


•  •  •  • 


106 


o  n  n 


LAMBDA=U*U/(32  2*(ZG-ZB)) 

BETA2=(8  0*U*U*(  1  0+DELTA*LAMBDA)**3 V(27  0*LAMBDA**2) 

IF  (BETA2  LT  0  0)  GO  TO  43 1 

BETA=SQRT(BETA2) 

XGB=U*U*BETA/32  2 
WRITE  (10,100)  U,XGB/LENGTH 
431  CONTINUE 
GO  TO  5000 

APPROXIMATE"  CUSP 

433  DO  434  1=1. JU 

WRITE  (*,2000)  I.JU 
U=UMIN+(UMAX-UMIN)*(1- 1  )/(JU- 1 ) 

DELTA=32  2*(ZW*MD-MW*ZD)/(ZD*WEIGHT) 

LAMBDA=U*U/(32  2*(ZG-ZB)) 

BETA2=(8  0*U*U*(1  0+DELTA*LAMBDA)**3)/27  0 
IF  (BETA2  LT  0  0)  GO  TO  434 
BETA=SQRT(BETA2) 

XGB=U*U*BETA/32  2 
WRITE  (10,100)  U.XGB/LENGTH 

434  CONTINUE 
C 

5000  STOP 

1001  FORMAT  ('  ENTER  1  SOLUTION  SETS  7 

1  '2  BIFURCATION  GRAPHS  ') 

1009  FORMAT  ('  ENTER  1  U  VARIATION  ’,/ 

1  '  2  XG  VARIATION ') 

1002  FORMAT  ('  ENTER  MINIMUM  VALUE  OF  U) 

1003  FORMAT  ('  ENTER  MAXIMUM  VALUE  OF  U’) 

1004  FORMAT  ('  ENTER  NUMBER  OF  INCREMENTS') 

1005  FORMAT  ('  ENTER  MINIMUM  VALUE  OF  XG/L') 

1006  FORMAT  ('  ENTER  MAXIMUM  VALUE  OF  XG/L') 

1010  FORMAT  ('  ENTER  VALUE  OF  XG/L') 

1012  FORMAT  ('  ENTER  VALUE  OF  U) 

1013  FORMAT  ('  ENTER  1  EXACT  COMPUTATION  ',/ 

1  '2  PITCHFORK  APPROXIMATION',/ 

2  '  3  SIMPLE  PITCHFORK  ',/ 

3  '  4  ANALYTIC  SET  ') 

1017  FORMAT  ('  ENTER  1  EXACT  CUSP  ',/ 

1  '2  APPROXIMATE  CUSP  ') 

2000  FORMAT  (215) 

100  FORMAT  (2E15  5) 


nno  onnn  n  ono 


END 


* 


C -  % 

SUBROUTINE  TRAP(NAB,OUT)  • 

NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE  ^ 

DIMENSION  A(1),B(1) 

NI=N-1  • 


OUT=0  0 
DO  1  I=LN1 

OUT  1  =0  5 *( A(I)+ A(I+ 1  ))*(B(I+ 1  )-B(I)) 
OUT=OUT+OUTl 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  EXS0LS(1VAR,L,IFC,U.XG,A1,A2,A3,A4,CDZ,AREA,ZW, 
&ZD.ZG,MW,MD,BUO,XA) 


IT  COMPUTES  AND  PRINTS  THE  EXACT  SOLUTION  SET  THETA 
(DEG  )  VERSUS  U  OR  XG 


REAL  MW.MD,PRNT,PI 
DIMENSION  VF(5,2) 

PI=4  0*ATAN(1  0) 

IF(IVAREQ  1)PRNT=U 
IF  (IVAR  EQ  2)  PRNT=XG 

FIND  FIRST  ESTIMATE  OF  SOLUTIONS 


•  • 


L=0 

VA=-89  5 
VA=VA*PI/180  0 
VAO=VA 

VO=THETEQ(  I  ,VA,  Al,  A2,A3,  A4,U) 

VAMIN=-89  5 
VAMAX=+89  5 
IVA=5000 
DO  10  I=2,IVA 
AI=I 

V  \=VAMIN+(  VAMAX- V  AMIN)*(I- 1  )/(IVA- 1 ) 
C  VA=AI-90  5 

VA=VA*PI/180  0 
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VAN=VA 

VN=THETEQ(  1 ,  VA,  A 1 ,  A2,  A3.  A4,U) 

VP=VO*VN 

IF  (VP  GE  O  0)  GO  TO  1 1 
L=L+1 

VF(L,l)=VAO 
VF(L,2)=VAN 
11  VO=VN 
VAO=VAN 
10  CONTINUE 

EXACT  COMPUTATION  OF  SOLUTIONS  VIA  NEWTON’S  METHOD 

E=1  E-8 
IEND=500 
DO  20  J=1,L 

X=(VF(J,l)+VF(J,2))/2  0 
C 

X1=X 

IXX=0 

IF  (IXX  EQ  0)  GO  TO  35 
C 

F=THETEQ(  1  ,X,  A 1 ,  A2.  A3,  A4,U) 
FDER=THETEQ(2,X,A1,A2,A3,A4,U) 

DO  30  K=1,IEND 
IF  (FDER  EQ  0  0)  STOP  1001 
DX=F/FDER 
X1=X-DX 

F=THETEQ(1,X1,A1,A2,A3,A4,U) 

FDER=THETEQ{2,X  1 ,  A1 ,  A2,  A3,  A4,U) 

IF  (F  EQ  0  0)  GO  TO  35 
A=ABS(X1-X) 

IF  (A-E)  35,35,40 
40  X=XI 
30  CONTITWE 
GOTO  20 
35  SOL=X1*180PI 
C  WRITE  (*,♦)  XI 
JPR=10+J 
WD=U*TAN(X1) 

DLT=((CDZ*AREA*WD*ABS(WD)HZW*U*WD))/(ZD*U*U) 

IF  ((ABS(DLT)  GT  0  4)  AND  (CDZ  NE  0  0))  THEN 


CALL 

SATUR(U.XG,CDZ,AREA,ZW,ZD,ZG,MW.MD.BUO,XA,DLT,PRNT) 

ELSEIF  ((ABS(DLT)  GT.O  4)  AND  (CDZ  EQ  0  0))  THEN 
CALL 

SATURl(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 

ENDIF 

IF  (L  EQ  1  AND  IFC  EQ  0)  WRITE  (10,100)  PRNT,SOL,DLT*  180/PI 
IF  (L  EQ  1 . AND  IFC.EQ  1 )  WRITE  ( 1 6, 100)  PRNT,SOL,DLT*  1 80/PI 
IF  (L  GT  1)  WRITE  (JPR,100)  PRNT.SOL.DLT*  180/PI 

20  CONTINUE 
RETURN 

100  FORMAT  (3E15, 5) 

END 

C - 

SUBROUTINE 

SATUR(U,XG,CDZ,AREA,ZW,ZD,ZG,MW,MD,BUO,XA,DLT,PRNT) 

REAL  MW,MD,PRNT 
C 

C  USED  TO  COMPUTE  SATURATION  CONDITIONS  FOR  THE  DIVE  PLANE 
WITH  DRAG 

C  COEFFICENT  INCLUDED 
C 

CD=CDZ 

IF  (DLT  GT  0  0)  THEN 

Wl=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U*0  4)) 
&  /(2*(-CD*AREA)) 

IF  (W1  GT  O.O)  THEN 

A=MW*U*W1-CD*XA*AREA*W1*W1-XG*BUO+MD*U*U*0  4 
B=.20*ZG*BUO 

C=MW*U*W1-CD*XA*AREA*W1*W1+XG*BUO+MD*U*U*0  4 
X 1  =(-B+SQRT((B  1  *B  1  )-4*(A*C)))/(2  0* A) 

SOL1=2  0*ATAN{X1) 

X2=(-B-SQRT((B  I  *B  1  )-4*(A*C)))/(2  0* A) 

SOL2=2.0*ATAN(X2) 

WRITE  (21,200)  PRNT,SOLl,SOL2 
ENDIF 

W2=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U*0  4)) 
&  /(2*(-CD*AREA)) 

IF  (W2  GT  0  0)  THEN 

A=MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO+MD*U*U*0  4 
B=-20*ZG*BUO 

C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*0  4 
X 1  =(-B+SQRT((B  1  *B  1  )-4*( A*C)))/(2  0* A) 


SOL1=2  0*ATAN(X1) 

X2=(-B-SQRT((B  1  *B  1  )-4*( A*C)))/(2  0* A) 

SOL2=2  0*ATAN(X2) 

WRITE  (22,200)  PRNT,SOLl.SOL2 
ENDIF 

W3=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0  4)) 
/(2*(CD*AREA)) 

IF(W3  LT  O  O)  THEN 

A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*0  4 
B=-20*ZG*BUO 

C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*0  4 
Xl=(-B+SQRT((B1*B1)-4*(A*C)))/(2  0*A) 

SOL1=2  0*ATAN(X1) 

X2=(-B-SQRT((B  1  *B  I  )-4*(  A*C)))/(2,0*  A) 

SOL2=2  0*ATAN(X2) 

WRITE  (23,200)  PRNT,SOLl,SOL2 
ENDIF 

W4=(-ZW*U-SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U*0  4)) 
/(2*(CD*AREA)) 

IF  (W4  LT  0  0)  THEN 

A=MW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*0  4 
B=-2  0*ZG*BUO 

C=MW*U*W4+CD*XA*AREA*W4*W4+XG*BUO+MD*U*U*0  4 
X 1  =(-B+SQRT((B  1  *B  1  )-4*( A*C)))/(2  0* A) 

SOL1=2  0*ATAN(X1) 

X2=(-B-SQRT((B  I  ‘B 1  )-4*( A*C)))/(2  0* A) 

SOL2=2  0*ATAN(X2) 

WRITE  (24,200)  PRNT,SOLl,SOL2 
ENDIF 

ELSE 

Wl=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(-CD*AREA)*ZD*U*U* 

(-0  4)))/(2*(-CD*AREA)) 

IF  (W1  GTO  O)  THEN 

A=MW*U*W1-CD*XA*AREA*WI*W1-XG*BUO+MD*U*U*(-0  4) 
B=-20*ZG*BUO 

C=MW*U*W1-CD*XA*AREA*W1*W1+XG*BUO+MD*U*U*(-0  4) 
X 1  =(-B-^SQRT((B  1  *B  1  )-4*( A*C)))/(2  0* A) 

SOL1=2  0*ATAN(XI) 

X2=(-B-SQRT((BI*Bl)-4*(A*C)))/(2  0*A) 

SOL2=2  0*ATAN(X2) 

WRITE  (3 1 ,200)  PRNT,SOL  1  .SOL2 
ENDIF 

W2=(-ZW*U-SQRT((ZW*U)*(ZWU)-4*(-CD*AREA)*ZD*U*U* 


&  (-0  4)))  /(2*(-CD*AREA)) 

IF  (W2  LT  0  0)  THEN 

A=MW*U*W2-CD*XA*AREA*W2*W2-XG*BUO+MD*U*U*(-0  4) 
B=-20*ZG*BUO 

C=MW*U*W2-CD*XA*AREA*W2*W2+XG*BUO+MD*U*U*(-0  4) 
X 1  =(-B+SQRT((B  1  *B  1  )-4*( A*C)))/(2  0*  A) 

SOL1=2  0*ATAN(X1) 

X2=(-B-SQRT((B  1  ‘B 1  )-4*( A*C)))/(2  0* A) 

SOL2=2  0*ATAN(X2) 

WRITE  (32,200)  PRNT,SOL21,SOL22 
ENDIF 

W3={-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U* 

&  (-0  4)))/(2*(CD*AREA)) 

IF  (W3  LT  O  O)  THEN 

A=MW*U*W3+CD*XA*AREA*W3*W3-XG*BUO+MD*U*U*(-0  4) 
B=-2  0*ZG*BUO 

C=MW*U*W3+CD*XA*AREA*W3*W3+XG*BUO+MD*U*U*(- 

0  4) 

X 1  =(-B+SORT((B  1  *B  1  )-4*(A*C)))/(2  0*  A) 

SOL1=2  0*ATAN(XI) 

X2=(-B-SQRT((B  I  *B  1  )-4*(A*C)))/(2  0* A) 

SOL2=2  0*ATAN(X2) 

WRITE  (33,200)  PRNT,SOLl,SOL2 
ENDIF 

W4=(-ZW*U+SQRT((ZW*U)*(ZW*U)-4*(CD*AREA)*ZD*U*U* 

&  (-0  4)))/(2*(CD*AREA)) 

IF  (W4  LT  0  0)  THEN 

A=NfW*U*W4+CD*XA*AREA*W4*W4-XG*BUO+MD*U*U*(-0  4) 
B=-2  0*ZG*BUO 

C=MW*U*W4+CD*XA*AREA*W4*W4+XG*BUO+MD*U*U*t- 

04) 

X 1  =(-B+SQRT((B I  *B  1  )-4*( A*C)))/(2  0* A) 

SOL  1=2  0*ATAN(X1) 

X2=(-B-S0RT((B  1  *B  I  )-4*( A*C)))/(2  0* A) 

SOL2=2  0*ATAN(X2) 

WRITE  (34.200)  PRNT.SOLI,SOL2 
ENDIF 
ENDIF 

200  FORMAT  (3E 15  5) 

END 

C . 

SUBROUTINE 

S  ATUR 1  (U.XG.CDZ.  AREA.ZW.ZD,ZG,MW,MD,BUO.XA,DLT,PRNT) 


o  n  o  o  n  n  n 


REAL  MW,MD,PRNT 


C 

C  COMPUTES  DIVE  PLANE  SATURATION  NEGLECTING  DRAG 
C 

PI=4  0*ATAN(I  0) 

A1=MD*U*U*0  4-MW*U*U*ZD*0  4-XG*ZW*BUO 
A2=-MD*U*U*0  4+MW*U*U*ZD*0  4-XG*ZW*BUO 
B1=-2  0*ZG*ZW*BUO 

C 1  =-MW*U*U*ZD*0  4+XG*ZW*BUO+MD*U*U*ZW*0  4 
C2=MW*U*U*ZD*0  4+XG*ZW*BUO-MD*U*U*ZW*0  4 
XCDl=((-Bl/Al)+SQRT(((Bl/Al)**2)-4  0*(Cl/Al)))/2  0 
XCD2=((-B  1  /A1  )-SQRT(((B  1/A  1  )**2)-4  0*(C  1  /A1  )))/2.0 
XCD3=((-Bl/A2)+SQRT(({Bl/A2)**2)-4  0*(C2/A2)))/2  0 
XCD4=((-BI/A2)-SQRT(((Bl/A2)**2)-4  0*(C2/A2)))/2.0 
XC 1  =2  0*ATAN(XCD1)*  180/PI 
XC2=2  0*ATAN(XCD2)*180/P1 
XC3=2  0*ATAN(XCD3)*  180/PI 
XC4=2  0*ATAN(XCD4)*  180/PI 
WRITE  (20,100)  PRNT,XC1,DLT*  180/PI 
WRITE  (41,100)  PRNT,XC2,DLT*  180/PI 
WRITE  (42,100)  PRNT,XC3,DLT*  180/PI 
WRITE  (43,100)  PRNT,XC4,DLT*  180/PI 
100  FORMAT  (3E15  5) 

END 

C . 

SUBROUTINE  PITSLS(IVAR,LL,IFC,U,XG,ALA2,A3,A4.A5) 

IT  COMPUTES  AND  PRINTS  THE  PITCHFORK  SOLUTION  SET  THETA 
(DEG  )  VERSUS  U  OR  XG 

DIMENSION  VF(5,2) 

PI=4  0*ATAN(1  0) 

IFdVAREQ  1)PRNT=U 
IF(IVAR  EQ  2)  PRNT=XG 

FIND  FIRST  ESTIMATE  OF  SOLUTIONS 


L=0 

VA=-89  5 
VA=VA*PI/I80  0 
VAO=VA 

VO=PlTCEQ(  I ,  VA,  A1 ,  A2,  A3,  A4,  A5,U) 
DO  10  1=2,180 


o  n  o 


AI=I 

VA=AI-90.5 

VA=VA*PI/180.0 

VAN=VA 

VN=PITCEQ(  1  ,VA,  A1 ,  A2,  A3,A4,A5,U) 

VP=VO*VN 

IF  (VP  GE  O  0)  GO  TO  11 
L=L+1 

VF(L,l)=VAO 
VF(L.2)=VAN 
11  VO=VN 
VAO=VAN 
10  CONTINUE 

EXACT  COMPUTATION  OF  SOLUTIONS  VTA  NEWTON’S  METHOD 


* 


«/ 


Ar 


E=l.E-5 
IEND=500 
DO  20  J=1,L 

X=(VF(J,l)+VF(J,2))/2  0 
F=PITCEQ(1,X,A1,A2,A3,A4,A5,U) 

X1=X 

IXX=0 

IF  (IXX  EQ  0)  GO  TO  35 

FDER=PITCEQ(2,X,A1,A2,A3,A4,A5,U) 

DO  30K=1,IEND 
IF  (FDER  EQ  0  0)  STOP  1001 
DX=F/FDER 
X1=X-DX 

F=PITCEQ(  1  ,X  1 ,  A1 ,  A2,  A3,  A4,  A5,U) 
FDER=PITCEQ(2,X  1 ,  A1 ,  A2,  A3,A4,A5,U) 

IF  (F  EQ  0  0)  GO  TO  35 
A=ABS(X1-X) 

IF  (A-E)  35,35,40 
40  X=X1 
30  CONTINUE 
GOTO  20 

35  SOL=Xl*  180  0/PI 
JPR=10+J 

IF  (L  EQ  1  AND  IFC  EQ  0)  WRITE  (10,100)  PRNT,SOL 
IF  (L  EQ  1  AND  IFC  EQ  1)  WRITE  (16,100)  PRNT,SOL 
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•  • 


non  ooo  noHOO 


WRITE  (JPR,100)  PRNT.SOL 


IF(LGT  1) 

20  CONTINUE 
RETURN 

100  FORMAT  (2E 15  5) 
END 


C - 

SUBROUTINE  SIMSLS(IVAR,L.IFC.U,XG.A1.A3.A4) 


IT  COMPUTES  AND  PRINTS  THE  SIMPLE  PITCHFORK  SOLUTION  SET 
■lETA 

(DEG  )  VERSUS  U  OR  XG 


DIMENSION  VF(5,2) 

PI=4  0*ATAN(1  0) 

IF  (IVAR  EQ  1)PRNT=U 
IF  (IVAR  EQ  2)  PRNT=XG 

FIND  FIRST  ESTIMATE  OF  SOLUTIONS 


L=0 

VA=-89,5 

VA=VA*PI/1800 

VAO=VA 

VO=SIMPEQ(  1 ,  VA  A1 ,  A3,  A4,U) 

DO  101=2,180 
AI=I 

VA=AI-90  5 
VA=VA*PI/180  0 
VAN=VA 

VN=SIMPEQ(1,VAA1,A3,A4,U) 

VP=VO*VN 

IF  (VP  GE  O  0)  GO  TO  1 1 
L=L+1 

VF(L.l)=VAO 
VT(L,2)=VAN 
11  VO=VN 
VAO=VAN 
10  CONTINUE 

EXACT  COMPUTATION  OF  SOLUTIONS  VIA  NEWTON'S  METHOD 

E=1  E-5 
IEND=500 


•  •  •  • 


•  •  •  •  • 
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DO20J=l.L 

X=(VF(J,l)+VF(J,2))/2.0 

C 

X1=X 

IXX=0 

IF  (IXX  EQ  0)  GO  TO  35 
C 

F=SIMPEQ(  1  ,X,  A 1 .  A3,  A4,U) 

FDER=SIMPEQ(2,X,A1,A3.A4.U) 

DO  30  K=  FIEND 
IF  (FDER  EQ  0  0)  STOP  1001 
DX=F/FDER 
X1=X-DX 

F=SIMPEQ(  1  ,X  1 ,  A 1 ,  A3,  A4.U) 

FDER=SIMPEQ(2,X  1 ,  A1 ,  A3,  A4.U) 

IF  (F  EQ  0  0)  GO  TO  35 
A=ABS(X1-X) 

IF  (A-E)  35,35,40 
40  X=X1 
30  CONTINUE 
GOTO  20 

35  SOL=Xl*  180  0/PI 
PR=10+J 

IF(L  EQ  1  AND  IFC  EQ  O)  WRITE  (10,100)  PRNT,SOL 
IF(L  EQ  1  AND  IFC  EQ  l)  WRITE  (16,100)  PRNT,SOL 
IF  (L  GT  1 )  WRITE  (JPR,  100)  PRNT,SOL 

20  CONTINUE 
RETURN 

100  FORMAT  (2E15  5) 

END 

C - 

SUBROUTINE  EXNUMS(NUM,A1,A2,A3,A4,U) 

C 

C  IT  COMPUTES  THE  NUMBER  OF  SOLUTIONS  OF  THE  EXACT  SOLUTION 

SET 

C 

PI=4  0*ATAN(l  0) 

L=0 

VA=-89  5 
VA=VA*PI/180  0 

V0=THETEQ(1,VA.A1,A2,A3,A4,U) 

DO  10  1=2,180 


onoon  no«oo 


VA=AI-90  5 
VA=AI*PI/180  0 

VN=THETEQ(  1  ,VA,  A1 ,  A2,A3,A4.U) 
VP=VO*VN 

IF  (VP  GE  O.  0)  GO  TO  1 1 
L=L+1 
11  VO=VN 
10  CONTINUE 
NUM=L 
RETURN 
END 


SUBROUTINE  PINUMS(NUM,A1,A2,A3.A4,A5,U) 


IT  COMPUTES  THE  NUMBER  OF  SOLUTIONS  OF  THE  PITCHFORK 
OLUTION  SET 


PI=4.0*ATAN(1.0) 

L=0 

VA=-89  5 
VA=VA*PI/1800 

V0=PITCEQ(1,VA,A1,A2,A3,A4,A5,U) 

DO  101=2,180 
AI=I 

VA=AI-90  5 
VA=VA*PI/180  0 

VN=PITCEQ(  1 ,  VA,  A1  ,A2,A3,A4,A5,U) 
VP=VO*VN 

IF  (VP  GE  0  0)  GO  TO  11 
L=L+1 
11  VO=VN 
10  CONTINUE 
NUM=L 
RETURN 
END 


SUBROUTINE  SINUMS(NUM,A1,A3,A4.U) 

IT  COMPUTES  THE  NUMBER  OF  SOLUTIONS  OF  THE  SIMPLE  PITCHFORK 
SOLUTION  SET 

PI=4  0*ATAN(1  0) 

L=0 


VA=-89  5 

VA=VA*Piyi80.0 

VO=SIMPEQ(  1 ,  VA,  A1 ,  A3,  A4,U) 

DO  101=2,180 
AI=I 

VA=AI-90.5 
VA=VA*PI/180.0 
VN=SIMPEQ(  1  ,VA,  A1  ,A3,  A4,U) 

VP=VO*VN 

IF  (VP.GE  0.0)  GO  TO  1 1 
L=L+1 
11  VO=VN 
10  CONTINUE 
NUM=L 
RETURN 
END 

C - 

FUNCTION  THETEQ(K,XA,A1,A2,A3,A4,U) 

C 

C  IT  COMPUTES  THE  VALUE  OF  THE  EXACT  EQUATION  FOR 
C  THETABAR  FOR  A  GIVEN  VALUE  OF  THETA 
C 

C  K  =  1  COMPUTE  THE  VALUE  OF  THE  FUNCTION 
C  K  =  2  COMPUTE  THE  VALL'E  OF  ITS  DERIVATIVE 
C 

SP=SIN(XA) 

CP=COS(XA) 

AP=ABS(SP) 

TP=TAN(XA) 

ATP=ABS(TP) 

GO  TO  (10,20),  K 

10THETEQ=A1*U*U*SP*CP+A2*U*U*SP*AP+A3*CP**3+A4*SP*CP**2 
C  10  THETEQ=AI*U*TP+A2*U*U*TP*ATP+A3*CP+A4*SP 
GOTO  50 

20  THETEQ=A1  *U*U*CP*CP-A1  *U*U*SP*SP-3  0*A3*CP*CP*SP+A4*CP**3 
&  -2  0*A4*SP*SP*CP+2  0*A2*U*U*CP*AP 

50  RETURN 
END 

C - 

FUNCTION  PITCEQ{K,XA,A1,A2,A3,A4,A5,U) 

C 

C  IT  COMPUTES  THE  VALUE  OF  THE  PITCHFORK  EQUATION  FOR 
C  THETABAR  FOR  A  GIVEN  VALUE  OF  THETA 


T- 


«7 


♦ 


•  • 
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•  • 


n  n  n  n 


* 

K  =  1  COMPUTE  THE  VALUE  OF  THE  FUNCTION  “v 

K  =  2  COMPUTE  THE  VALUE  OF  ITS  DERIVATIVE  • 

TP  =TAN(XA)  ^ 

CP  =COS(XA) 

CP2=CP*CP 

AP  =ABS(TP)  • 

GO  TO  (10,20),  K 


10PITCEQ=A1*(U*TP)**3+A2*U*U*TP*AP+A3*U*TP+A4+A5*U*U*TP*TP 
GOTO  50 

20  PITCEQ=3  0*U*U*U*TP/CP2+2  0*A2*U*U*AP/CP2+A3*U/CP2 

&  +2  0*A5*U*U*TP/CP2  • 

50  RETURN 
END 

C - 

FUNCTION  SIMPEQ(K,XA,A1,A3,A4,U) 

C  • 

C  IT  COMPUTES  THE  VALUE  OF  THE  SIMPLE  PITCHFORK  EQUATION  FOR 
C  THETABAR  FOR  A  GIVEN  VALUE  OF  THETA 
C 

C  K  =  1  COMPUTE  THE  VALUE  OF  THE  FUNCTION 

C  K  =  2  COMPUTE  THE  VALUE  OF  ITS  DERIVATIVE  •  # 

C 

TP  =TAN(XA) 

CP  =COS(XA) 

CP2=CP*CP 

AP  =ABS(TP)  • 

GO  TO  (10,20),  K 

10  SIMPEQ=A1*(U*TP)**3+A3*U*TP+A4 
GOTO  50 

20  SIMPEQ=3  0*U*U*U*TP/CP2+A3*U/CP2 

50  RETURN  ^ 

END 
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%  THIS  IS  PROGRAM  BIFSU  M  IT  COMPUTES  THE  BIASED  BIFURCATION 
CUSP  FOR 

%  DELTA  SATURATED  INCLUDING  DRAG  TERMS  IT  USES  THE  FZERO 
FUNCTION  TO 
%  SOLVE  THE  EQUATIONS 

global  W  B 1  U  n  Zw  Zdlt  A  XA  CD  rho, 

%  ESTABLISH  GEOMETRIC  PARAMETERS 
W=1556.2363;B1=1.0001*W; 

A=19  8473;XA=0  20126; 

CD=0  3, 

Iy=561  32; 
g=32.2, 
m=W/g; 
rho=l  94, 

L=13  9792, 
xb=0, 

Alphas=l, 

Alphab=0, 
zg=0  1, 
zb=0, 

zgb=zg-zb, 

%  NON-DIMENSIONING  FACTORS 
ndl=  5*rho*L''2, 
nd2=  5*rho*L"3, 
nd3=  5*rho*L^4, 
nd4=  5*rho*L^5, 

%  HYDRODYNAMIC  COEFFICIENTS 

Zqdnd=-6  33e-4,Zwdnd=-l  4529e-2,Zqnd=7  545e-3,Zwnd=-I  39le-2, 

Zds=-5  603e-3.Zdb=0  5*(-5  603e-3),Zdltnd=(Alphas*Zds+AJphab*Zdb), 

Mqdnd=-8  8e-4,Mwdnd=-5  61e-4,Mqnd=-3  702e-3,Mwnd=l  0324e-2, 

Mds=-0  002409.Mdb=0  5*(0  002409),Mdltnd=(Alphas*Mds+Alphab*Mdb). 

Zqd=nd3*Zqdnd,Zwd=nd2*Zwdnd,Zq=nd2*Zqnd,Zw=ndl*Zwnd. 

Zdlt=ndl*Zdltnd. 

Mqd=nd4*Mqdnd,Mwd=nd3*Mwdnd,Mq=nd3*Mqnd,Mw=nd2*Mwnd, 

Mdlt=nd2*Mdltnd, 


%  ESTABLISH  SPEED  RANGE 
U=9  001  2  6, 
for  n=l  :length(U), 


%  SOLVES  FORCE  EQUATION  USING  FZERO  FUNCTION 
th(n)=fzero('satur',(-Zdlt/Zw),0.0000 1 ), 

%  COMPUTES  RELATIONSHIP  BETWEEN  U  AND  Xgb  WITH  RESULTS  FROM 
ABOVE  EQUATION 

xg(n)=(Mw*U(n)'^2*tan(th(n))-zg*W*sin(th(n))-(0  5)*rho*CD*A*XA*U(n)^2* 

tan(th(n))*abs(tan(th(n)))+Mdlt*U(n)''2*(0  4))/(W*cos(th(n))), 

end, 

XgbL=xg/13  9792, 

Lp=sqrt(U  ^2/(3  22)), 
for  n=l  length(U); 

th  1  (n)=fzero('satur  1  ',(-Zdlt/Zw),0  0000 1 ), 

xgl(n)=(Mw*U(n)^2*tan(thl(n))-zg*W*sin(thl(n))-(0  5)*rho*CD*A* 
XA*U(n)^2*tan(th  1  (n))*abs(tan(th  1  (n)))+Mdlt*U(n)^2*(-0  4))/(W*cos(th  1  (n))), 
end, 

XgbLl=xgl/13  9792, 


%  THIS  IS  PROGRAM  SATUR  M  IT  ESTABLISHES  THE  FUNCTION  TO  BE 
SOLVED  IN 

%  THE  BIFSU  M  PROGRAM 


function  y=satur(th), 

y=Zw*U(n)^2*tan(th)+(W-Bl)*cos(th)+Zdlt*U(n)''2*(0  4) 

-(0  5)*rho*CD*A*U(n)^2*tan(th)*abs(tan(th)),  * 

%  THIS  IS  PROGRAM  SATUR  I  M  IT  ESTABLISHES  THE  FUNCTION  TO  BE 
SOLVED  IN 

%  THE  BIFSU  M  PROGRAM 
function  yl=saturl(thl ), 

y  1  =Zw* U(n)" 2*tan(th  1  )-t-(W-B  1  )*cos(th  1  )+Zdlt *U(n)''2*(-0  4) 

-(0  5)*rho*CD*A*U(n)^2*tan(thl)*abs(tan(thI )). 
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APPENDIX  B 


DTRC  SUBOFF  MODEL  CHARACTERISTICS 


0BSI(2f  PASAMZTZRS  -  aUIi. 


;  L«n*th 

3«:vc«n  ?sr3«dlcular3 

(IC?) 

(f:)i 

13.9*9:  i 

lanjta 

Overall 

(LOA) 

(f;)l 

:.«ns:.n 

:3  :he  Cancer  of  Buoyancy 

C--C3) 

(f:)i 

6.60^2 

l«nK:a 

3f  ?or«oody 

t  •  \ 

(fOI 

3.3323 

of  Parallel  Jliddleboay 

(  f  *)l 

^  31*'f 

of  Run 

(f:)l 

3 . 

Oiaa«t«r 

(::>! 

1 . 566* 

rincness 

3.575 

Icefl 

Bare 

B.3.  - 

Mjjjjl 

3.3.  • 

•ully 

Bull 

Sail 

kSSQI 

Ring  Ving  1 

.Appended 

70Uf:^ 

24.692 

24.33434 

24.78279 

24.75678 

24.98991 

LCa(f:) 

6.59003 

6.5726 

6.612732 

6.50889 

•  5.513906 

7a(ft) 

0.0 

43.306711 

0.0 

0.3 

-0.306669 

7S(f:-) 

63.717 

65.354 

65.514 

63.71* 

57 . 551 

Buoyancy  (lb) 

1537 . 6341 

1546.5433 

1541.3380 

1541.-183 

1556.2262 

S' 

0.018078 

0.313182 

0.013144 

0.313122 

3.313296 

0.388697 

0.410569 

-0.111173 

-0.061225 

.  12**  *67 

-y' 

0.001053 

3.301059 

0.001066 

3.301066 

3.301C8<* 

Icca 

i  Fully 

1  Appended 

70L(f:^) 

2-. 9899 

lC3(f;) 

!  6.6139 

7a(f:) 

1  -0.006669 

n 

VS(f:-) 

1  67.651 

Buoyancy (lb) 

!  1556.2363 

:  0.313296 

t^'a'^lO* 

-0. 127467 

‘z' 

1  0.001084 

«7 


Ar 


^ondiMasianal  and  ztoss  aeccional  arsas  far  :.ie  lul’. . 


3Tat:qn 

!  3/3^ 

A/  A^ 

.V 

1 

:  0.00000 

0.00000 

0.1 

0.29033 

0 . 

0.2 

0.39396 

0.15520  ■ 

1  0.3 

0.46600 

0.21715 

I  0.- 

0.52147 

0.27194  ' 

1  0.3 

0.56627 

0.32066  , 

i  0.5 

0.50332 

0.36424 

i  0.7 

0.63314 

0.40340 

! 

0.70744 

0.50047  1 

1  2.0 

0.34713 
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